%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           
%       cake2.m: Cake eating problem II
%
%       The program solves the "Cake eating problem II" using the method
%       of undetermined coefficients.
%
%    NOTE:
%   1) To evaluate the residual you have to give a guess of the coefficients
%   2) To derive the optimal value funcion call the function with the minimization 
%
%     i.e. [coefficients,residual,conv,options]= fminsearch('dynport',[A,B])
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res=cake2(x)

% PARAMETERS
delta   = 0.5;                            % Exogenous return from savings
beta    = 0.5;                            % Discount rate

% Coefficients for the functional form for the value function v(k)= A + B log(k)
A    = x(1);
B    = x(2);

%----------------------------------------------------------------
% INITIALIZE VARIABLES AND FUNCTIONS
%----------------------------------------------------------------
% Define grid points for cake
nump  = 25;                             % Number of grid points
whigh = 5;                              % Highest point in the grid
wlow  = whigh/nump;                     % Lowest point in the grid  
w    = linspace(wlow,whigh,nump)';     % State space
tol   = 1e-5;                           % Tolerance criterium

% Guess for the initial value function and policy function
v1      = zeros(nump,1);                % Initialize value function
wprime  = zeros(nump,1);                % Initialize savings policy function
c    = zeros(nump,1);                % Initialize consumption policy function


%----------------------------------------------------------------
% COMPUTE VALUE FUNCTION
%----------------------------------------------------------------
% For every point in the grid solve the FOC and compute consumption
%  -U'(delta*w-w')+beta*V'(w')

    wprime = (delta*beta*B/(1+beta*B))*w;

    c = delta*w-wprime;

% Compute the value function associated to the optimal policy

    v1 = log(c) + beta* (A + B*log(wprime));

% Compute the value function associated to the guessed coefficients

    v0  = A + B*log(w);

%----------------------------------------------------------------
% COMPUTE THE RESIDUAL FUNCTION AND PLOT CONVERGENCE
%----------------------------------------------------------------

    res = max(abs(v1-v0));

% Plot convergence of the value function

	plot(w,v1,':rs',w,v0); 
	legend('v1','v0')   
	axis([wlow whigh -10 3])
    title('Plot of 2.c')
	pause(.1)