%%%%%%%%%%%%%%%%%%%%  Solow Model (1956)  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
% Steady State.                                                           %
%                                                                         %
% [1] Relations of Solow Model (1956):                                    %
%   dk/dt = s*f(k(t))-theta*k(t)                                          %
%   netVest = s*f(k)-theta*k                                              %
%   k(0) = k_0 > 0                                                        %
%                                                                         %
% [2] COBB-DOUGLAS FUNCTION:                                              %
%   Y = A*K^(alpha)*L^(beta)                                              %
%   with alpha>0, beta<1.                                                 %
%   If beta = 1 - alpha, constant returns to scale. Then the (per capita) %
%   cobb-Douglas function is given by:                                    %
%                      y = A*k^(alpha)                                    % 
%                                                                         %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                                                                         %
%                                                                         %
% see "SwSLID" "cobbdouglas.m", "netVestCB.m", "TnetVestCB.m", "VestCB"   %
%                                                                         %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global A alpha theta s;

% s = Savings rate
s=str2num(VAsc);

 A=10; 
 alpha=0.32; 
 theta=0.09; % Depreciation rate  

% Find k* for which netVest.m is zero
kstar=fzero(@netVestCB,[1 100],0.0001);

% Savings, Investment.
k1=linspace(0,100,100);
Inv=VestCB(k1);
dgn=theta*k1;

subplot(2,2,1);
xlabel(['Capital k. k* = ', num2str(kstar)]);
ylabel('Invest.');
plot(k1,Inv,'b',k1,dgn,'r');
xlabel(['Capital k. k* = ', num2str(kstar)]);
ylabel('Invest.');
title(['Solow. sf(k), \thetak, for s=', num2str(s)]);
grid;
hold on;
Dstar=theta*kstar;
plot(kstar,Dstar,'r.','markersize',20);
hold on;
lineStar=linspace(0,Dstar,30);
Qaxis=length(lineStar);
lineQ=ones(1,Qaxis).*kstar;
plot(lineQ,lineStar,'k-.');
hold off;

% Steady state
k=linspace(0,kstar,100);
NInv=netVestCB(k);
Nstar=netVestCB(kstar);

subplot(2,2,2);
plot(k,NInv); 
xlabel(['Capital k. k* = ', num2str(kstar)]);
ylabel('Net Invest.');
title('"sf(k)-\thetak" with a Cobb-Douglas');
grid;
hold on;
plot(kstar,Nstar,'r.','markersize',20);
hold off;

 Ini_Stat=1;
 Tini=0;
 Tf=100;

[t,y]=ode45(@TnetVestCB,[Tini Tf],Ini_Stat);

subplot(2,2,3);
plot(t,y);
xlabel('Time');
ylabel('Capital  (k_0=1)');
title('Time Evolution of sf(k(t))-\thetak(t)');
grid;
hold on;
edy=y(end,:);
plot(Tf,edy,'k.','markersize',20);
hold off;


 Ini_StatA =0.5;
 Ini_StatB =100;
 Tini=0;
 Tf=100; 
[tA,yA]=ode45(@TnetVestCB,[Tini Tf],Ini_StatA);
[tB,yB]=ode45(@TnetVestCB,[Tini Tf],Ini_StatB);

subplot(2,2,4);
plot(tA,yA);
hold on;
plot(0,Ini_StatA,'b.','markersize',20);
hold on;
plot(0,Ini_StatB,'b.','markersize',20);
hold on;
plot(tB,yB);
xlabel('Time');
ylabel('Capital');
title('Time Evolution for initial values k_0=1 k_0=100');
grid;
hold on;
edyA=yA(end,:);
edyB=yB(end,:);
plot(Tf,edyA,'k.','markersize',20);
hold on;
plot(Tf,edyB,'k.','markersize',20);
hold on;
lhA=length(tA);
lineA=ones(1,lhA).*kstar;
plot(tA,lineA,'k-.');
hold on;
lhB=length(tB);
lineB=ones(1,lhB).*kstar;
plot(tB,lineB,'k-.');
hold off


% Last revised on 9 March 2010.

%"Complex and Chaotic Nonlinear Dynamics. 
% Advances in Economics and Finance, 
% Mathematics and Statistics" 
% T.Vialar, Springer 2009.
% Copyright(c).