function F = steadystate_jh2(x)
global U1;global U2;global U3;global U4;global U5; global kout;global p;

n =48;
w = U1;
d0=x(1);
d1=x(2);
d2=x(3);
d3=x(4);


for i=1:n
p(i,1) = exp(-exp(d0 + d1*i + d2*i^2 + d3*i^3));
end

uhat(1,1) = w;

for j=2:n
uhat(j,1) = w*prod(p(1:j-1,1),1);
end
%[p uhat]
uhat1 = uhat(1,1);
uhat23 = sum(uhat(2:3,1),1);
uhat46 = sum(uhat(4:6,1),1);
uhat712 = sum(uhat(7:12,1),1);
uhat13 = sum(uhat(13:n,1),1);

X = [uhat1-U1;
     uhat23-U2;
     uhat46-U3;
     uhat712-U4;
     uhat13-U5];
if kout == 1
    [p uhat]
    X'
end
F = sqrt(X'*X);
F = F*1.e4;
