function x=gradhp(param);
a4=param(9); 
b4=param(10);
delta=param(13);
beta=param(14);
grate=param(15);
c0=param(16:20);
c1=zeros(5); c2=c1; c3=c1;
c1(:)=param(21:45); 
c2(:)=param(46:70);

x(6:10,1) = (eye(5)-c1-c2)\c0;
khm = ( (-a4 + ( (grate-beta*(1-delta)-beta*delta*x(9))/ ...
        (beta*a4*(1-x(9))) )^(b4/(1-b4)) )/(exp(x(7))^b4*(1-a4)) )^(-1/b4);
z0=[250;300;25];
z = newton('sshp4b',z0,[param(1:15);x(6:10);khm],1e-7,100);

hm=z(1); hn=z(2); khn=z(3);
x(1) = khm*hm;
x(2) = hm;
x(4) = hn;
x(5) = khn*hn;
x(3) = (grate-1+delta)*(x(1)+x(5));
x(11:12) = x(1:2);
