function Q_Flow = PressFlowRelation(xu,xo,P2,P3,Qin,par)
x0 = [P2,P3];
if xo<0.1
    xo=0.05;
end

Cv1=par.Aup*par.Cdu;%Pipe dia underflow
Cv2=par.Aop*par.Cdo;%pipe dia overflow

%% Qin as a sum of Qo and Qu
%options = optimoptions('fsolve','Display','off');

options = optimoptions('fsolve','Display','off','OptimalityTolerance',1e-10,'FiniteDifferenceType','central');
x = fsolve(@(x) PressEquations_withQin(Qin,xu,xo,par,x),x0,options);
Qu=abs(Cv1*xu*sqrt(2/par.Rho_w*(x(2)-par.Patm)));
Qo=abs(Cv2*xo*sqrt(2/par.Rho_o*(x(1)-par.Patm)));
P2=abs(x(1));
P3=abs(x(2));
Q_Flow =[Qu Qo P2 P3];