%%%%%%%%%%%%%%%%%%%%%%%%%
%% steady state solver %%
%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

% Average (1976m01-2013m12):  U1=3.210;U2=2.303;U3=1.238;U4=1.050;U5=0.636;
% Average (1976m01-2017m06):  U1=3.178;U2=2.281;U3=1.244;U4=1.064;U5=0.664;
% Average (2007m12-2013m12):  U1=3.339;U2=2.787;U3=2.131;U4=2.426;U5=1.902; 

Udata(1,:) = [3.178,2.281,1.244,1.064,0.664];
Udata(2,:) = [3.339,2.787,2.131,2.426,1.902];
 
for i=1:2
    U1 = Udata(i,1);    U2 = Udata(i,2);    U3 = Udata(i,3);
    U4 = Udata(i,4);    U5 = Udata(i,5);
    
%% 5eqn/5unkown %%

global U1;global U2;global U3;global U4;global U5; global kout;global p;
x0 =[0.1;0.1;0.1;0.1];
lb =[-1;-1;-1;-1];
ub =[1;1;1;1];
kout = 0;
%gs = GlobalSearch;
%problem = createOptimProblem('fmincon','x0',x0,...
%    'objective',@steadystate4,'Aineq',[],'bineq',[],'lb',lb,'ub',ub);
%[x,fval,flag] = run(gs,problem);

[x,fval,flag] = fmincon(@steadystate_jh2,x0,[],[],[],[],lb,ub);
kout = 1;
steadystate_jh2(x);
kout = 0;

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

nn = 2*i;
table_jh(nn,:) = [w,d0,d1,d2,d3,U1,U2,U3,U4,U5];
P(:,i)=p;
end
format long;
table_jh
format short;