function f = drillsolve4u(y)
%Function that solves system of drilling equations.
%Run main4u.m
%DEH - 2010

global db Aa dp g Wdbt R0 rho rhos theta1 D WOB RPM qmud pc

f = zeros(3,1);

newrho = rho*(1-y(2)/100) + rhos*y(2)/100; %rhof, kg/m3
bhp = pc + theta1*qmud/1000 + newrho*g*D/1e5; %BHP, bar

An = 3*pi*0.01^2/4; %nozzle x-sec area, m2
Fj = qmud^2*rho/An/1e6; %hydraulic jet impact force, N
ROP = R0*(WOB/db-Wdbt)/(71.433-Wdbt)*(RPM/60)^0.7*exp(0.01*...
    (470-bhp))*(Fj/4448.22)^.3;

va = qmud/1000/Aa; %annulus velocity, m/s
vslip = sqrt(8/9*g*dp*(rhos-newrho)/rho); %slip velocity, m/s
vT = va - vslip; %transport velocity, m/s
qs = ROP*pi/4*db^2*1000/3600; %cuttings feed rate, l/s
xc = qs/1000/Aa/vT*100; %cuttings fraction

f(1) = ROP - y(1);
f(2) = xc - y(2);
f(3) = bhp - y(3);