global m n nr co Ir Hr cL
global rho mu grav
m=8;n=7;nr=2;
co=[1 2;2 4;4 3;4 6;4 5;3 5;6 5;6 7];
Ir=[1 7];Hr=[200 210];
cL=-2*[0 0 0.01 0.01 0.01 0.01 0];
rho=1000;mu=1e-3;grav=9.8;
q=0.01*ones(m,1);
h=[200*ones(n-1,1);210];
vR=0*ones(nr,1);
x0=[q;h;vR];
res=fhidr1(x0)
[x fval info]=fsolve(@fhidr1,x0)