%Inintial values
BETA=0.97;
DELTA=0.05;
ALFA=0.3;
THETA=1.5;
DOWN=0.6;



syms  b d e f g h l m s
eq1=b-DELTA;
eq2=d-(1+THETA);
eq3=e-1/BETA;
eq4=f-(1-DELTA);
eq5=g-ALFA;
eq6=h-THETA;
eq7=l-1/(ALFA-1);
eq8=m-ALFA/(ALFA-1);
eq9=s-DOWN;
[b d e f g h l m s]=solve(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9);


b1=eval(b(1));
d1=eval(d(1));
e1=eval(e(1));
f1=eval(f(1));
g1=eval(g(1));
h1=eval(h(1));
l1=eval(l(1));
m1=eval(m(1));
s1=eval(s(1));


syms w
eq=b1*d1*(w^d1)*(((e1-f1)*(w^d1)/g1/(h1*((w^d1)/d1)+1/d1))^l1)/h1/(w^d1-1)-g1*(((e1-f1)*(w^d1)/g1/(h1*((w^d1)/d1)+1/d1))^m1)-s1;
[w]=solve(eq);
j=find(abs(eval(w))>1);
x=eval(w(j));
k=((1/BETA-(1-DELTA))*(x^(1+THETA))/ALFA/(THETA*((x^(1+THETA))/(1+THETA))+1/(1+THETA)))^(1/(ALFA-1));
y=k^ALFA;
i=(ALFA*y+DOWN)*(1-x^(-THETA));