% Define parameters
syms alp delt bet sig mu ra  rg  as gs gsy  ys
%
% Define variables 
syms k1 a1 g1 c1 lam1 w1 xx1 yy1 h1 r1 l1
syms k a g c lam w xx yy h r l 


alp     = 0.35;
delt	= 0.025;
bet     = 0.99;
sig     = 1.500;
mu      = 0.300;
gsy		= 0.200;
ra		= 0.9500;
rg		= 0.97;
as      =   1;

% Function definitions
ut=(c^mu*(1-h)^(1-mu))^(1-sig)/(1-sig);
pr=a*k^alp*h^(1-alp);     
%
e1	= k1-xx-(1-delt)*k;
%
e2	=  log(a1)-(1-ra)*log(as)-ra*log(a);
%
e3	= log(g1)-(1-rg)*log(gs)-rg*log(g);
%
e4	= diff(ut,'c')-lam;
%
e5 = diff(ut,'h')+lam*w;
%
e6	= lam-bet*lam1*(r1+1-delt);
%
e7	= w-diff(pr,'h'); 
%
e8	= r-diff(pr,'k'); 
%
e9	= yy-pr;
%
e10	= yy-c-xx-g;
%
f=[e1;e2;e3;e4;e5;e6;e7;e8;e9;e10];
%
x = [k a g] ;
y = [c h lam w r yy xx];
x1 = [k1 a1 g1] ;
y1 = [c1 h1 lam1 w1 r1 yy1 xx1];
nx = length(x);
ny = length(y);

%Make f a function of the logarithm of the state and control vector

f = subs(f, [x,y,x1,y1], exp([x,y,x1,y1]));

%Compute analytical derivatives of f

fx = jacobian(f,x);
fx1 = jacobian(f,x1);
fy = jacobian(f,y);
fy1 = jacobian(f,y1);

%Get steady state values

% Initial guesses for the steady state values
 
k0=1.5; 				%capital
c0=.8;				%consumption
h0=0.3; 				%work
xx0=.1; 				%investment
r0=0.04;				%interest rate
w0=2; 				%wage
y0=1.5;
lam0=0.2;

x0=[k0;c0;h0;xx0;r0;w0;y0;lam0]; 	
sol=fcsolve('rbc_steady',x0,[])			          % Call the function 'rbc_steady'
													% which solves the system of steady state 
													% equations. The solutions are in the vector sol 
% Steady state values: Assignment

ys=sol(1);  
cs=sol(2);  
ks=sol(3);  
hs=sol(4); 
xxs=sol(5);  
rs=sol(6); 
ws=sol(7);  
lams=sol(8);

gs = gsy*ys;
as = 1;


c  = log(cs);
c1  = log(cs);
xx  = log(xxs);
xx1  = log(xxs);
yy  = log(ys);
yy1  = log(ys);
h  = log(hs);
h1  = log(hs);
k  = log(ks);
k1  = log(ks);
w  = log(ws);
w1  = log(ws);
r  = log(rs);
r1  = log(rs);
lam  = log(lams);
lam1  = log(lams);
a = log(as);
a1 = log(as);
g = log(gs);
g1 = log(gs);

nfx = zeros(size(fx));
nfx(:) = eval(fx(:));

nfx1 = zeros(size(fx1));
nfx1(:)= eval(fx1(:));

nfy = zeros(size(fy));
nfy(:) = eval(fy(:));

nfy1 = zeros(size(fy1));
nfy1(:)= eval(fy1(:));

nf = zeros(size(f));
nf(:)=eval(f(:));


A = [-nfx1 -nfy1];

B = [nfx nfy];

[cx,sx]=solab(A,B,size(nfx,2));