% Define parameters
syms alp delt bet gam rho ks zs cs
%
% Define variables 
syms k1 z1 c1
syms k z c 

alp	= 0.35;
delt =0.025;
gam	= 2;
bet	= 0.988;
rho	= 0.9;


% Compute the steady state manually 
zs = 1 ; %steady-state value of technology shock 

ks = ((1/bet+delt-1)/alp)^(1/(alp-1));    %steady-state value of capital
cs = zs*ks^(alp)-delt*ks;                 %steady-state value of consumption 


e1 = c + k1- (1-delt)*k - z*k^alp;
e2 = log(z1) - rho*log(z);
e3 = c^(-gam) - bet*c1^(-gam)*(z1*alp*k1^(alp-1)+1-delt);
%
f=[e1;e2;e3];
%
x = [k z] ;
y = [c ];
x1 = [k1 z1] ;
y1 = [c1];

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
 
k=log(ks); 				%capital
k1=log(ks); 				%capital
c=log(cs);				%consumption
c1=log(cs);
z=log(zs); 				%techno
z1=log(zs);

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));