%HJB EQUATION FOR NEOCLASSICAL GROWTH MODEL clear all; clc; tic; s = 2; a = 0.3; d = 0.05; r = 0.05; A = 1; kss = (a*A/(r+d))^(1/(1-a)); css = A.*kss.^a - d.*kss; dVss = css.^(-s); I=151; kmin = 0.5*kss; kmax = 1.5*kss; k = linspace(kmin,kmax,I)'; dk = (kmax-kmin)/(I-1); iss = (I+1)/2; check = k(iss)-kss; %CHECK THAT MID-POINT OF GRID IS kss disp(check) Delta = 0.2*dk; %CANNOT BE TOO LARGE RELATIVE TO dk, SEE CANDLER p.181 (DISCUSSION OF CFL NUMBER) maxit=10000; crit = 10^(-6); dV = zeros(I,1); c = zeros(I,1); %INITIAL GUESS v(:,1)= (A.*k.^a - d.*k).^(1-s)/(1-s)/r; for n=1:maxit V = v(:,n); %AT STEADY STATE: USE BOUNDARY CONDITION dV(iss) = dVss; %BELOW STEADY STATE: USE FORWARD DIFFERENCE for i=1:iss-1 dV(i)=(V(i+1)-V(i))/dk; end %ABOVE STEADY STATE: USE BACKWARD DIFFERENCE for i=iss+1:I dV(i)=(V(i)-V(i-1))/dk; end c = dV.^(-1/s); %FOC U = c.^(1-s)/(1-s); Vchange(:,n) = U + dV.*(A.*k.^a - d.*k - c) - r.*V; v(:,n+1) = v(:,n) + Delta*Vchange(:,n); dist(n) = max(abs(Vchange(:,n))); if dist(n)