function f = drillsolve(y)
%Function that solves system of drilling model equations.
%Run script main.m.
%DEH - 2010

global db Aa dp g Wdbt R0 rho rhos WOB RPM qmud bhp

f = zeros(2,1);

newrho = rho*(1-y(2)/100) + rhos*y(2)/100;

An = 3*pi*0.01^2/4; %Nozzle x-sec area, m^2
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 = y(1)*pi/4*db^2*1000/3600; %Cuttings feed rate, l/s
xc = qs/1000/Aa/vT*100; %Cuttings percent

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