function [f,fj]=sshp4b(z,p) 

H=1134;
hm=z(1);  hn=z(2);    khn=z(3);

a0=p(1);  b0=p(2);    a1=p(3);   b1=p(4);  a2=p(5); b2=p(6);
a3=p(7);  b3=p(8);    a4=p(9);   b4=p(10); b=p(11);
r =p(12); delta=p(13);  beta=p(14); mu=p(15); g=p(16); 
sm=p(17); sn=p(18); tk=p(19); tn=p(20); khm=p(21); 
del=mu-1+delta;
a41 = (1-a4)*exp(sm)^b4;
a31 = (1-a3)*exp(sn)^b3;

i    = del*(khm*hm+khn*hn);
cm   = (1-g)*(a4*khm^b4+a41)^(1/b4)*hm-i;
cn   = (a3*khn^b3+a31)^(1/b3)*hn;
h    = a0*hm^b0+(1-a0)*hn^b0; 
cg   = g*(a4*khm^b4+a41)^(1/b4)*hm;
cp   = a2*cm^b2+(1-a2)*cn^b2;
c    = a1*cp^(b1/b2)+(1-a1)*cg^b1;       
L    = H-2*h^(1/b0);

dcm  = [(cm+i)/hm-del*khm;-del*khn;-del*hn];
dcn  = [0;cn/hn;a3*khn^(b3-1)*cn/(a3*khn^b3+a31)];
dcp  = a2*b2*cm^(b2-1)*dcm+(1-a2)*b2*cn^(b2-1)*dcn;
dhm  = [1;0;0];
dhn  = [0;1;0];
dkhn = [0;0;1];
dcg  = cg*dhm/hm;
dc   = a1*b1/b2*cp^(b1/b2-1)*dcp+(1-a1)*b1*cg^(b1-1)*dcg;
dh   = a0*b0*hm^(b0-1)*dhm+(1-a0)*b0*hn^(b0-1)*dhn;
dL   = -2*h^(1/b0-1)/b0*dh;

f1 = c/L*hm^(b0-1)*h^(1/b0-1)*cm^(1-b2)*cp^(1-b1/b2);
f(1,1)=f1-a1*a2*a41*(1-tn)*b*(a4*khm^b4+a41)^(1/b4-1)/(2*a0*(1-b));

f2 = c/L*hn^(b0-1)*h^(1/b0-1)*cn^(1-b2)*cp^(1-b1/b2)* ...
     (a3*khn^b3+a31)^(1-1/b3);
f(2,1)=f2-a1*(1-a2)*a31*b/(2*(1-a0)*(1-b));
 
f3 = khn^(b3-1)*(a3*khn^b3+a31)^(1/b3-1)*(cm/cn)^(1-b2);
f(3,1)= f3-a2*(mu-beta+beta*delta)/(beta*a3*(1-a2));

%Y=cg/g
%str1='k,hm,i,y,g,kn,taxn'
%ss=[khm*hm+khn*hn,hm,i,Y,g*Y,khn*hn,tn] 
%str2='km,hn,cn/y,cm/y,kn/(km+kn),km/y'
%ss=[khm*hm,hn,cn/Y,cm/Y,khn*hn/(khn*hn+khm*hm),khm*hm/Y]

if nargout > 1; 
  fj(1,:) = f1*(dc/c-dL/L+(b0-1)*dhm/hm+(1/b0-1)*dh/h+(1-b2)*dcm/cm+ ...
            (1-b1/b2)*dcp/cp)';
  fj(2,:) = f2*(dc/c-dL/L+(b0-1)*dhn/hn+(1/b0-1)*dh/h+(1-b2)*dcn/cn+ ...
            (1-b1/b2)*dcp/cp)';
  fj(2,3) = fj(2,3) + f2*(b3-1)*a3*khn^(b3-1)/(a3*khn^b3+a31);
  fj(3,:) = f3*(((b3-1)/khn+(1-b3)*a3*khn^(b3-1)/(a3*khn^b3+a31))*dkhn+ ...
            (1-b2)*(dcm/cm-dcn/cn))';

end;

