function solve_two_period_problem()

% Assigning parameters
crra            = 0.5  ;
discount_factor = 0.99 ;
interest_rate   = 1.04 ;
initial_wealth  = 1    ;

% Packaging parameters
parameters.crra            = crra            ;
parameters.discount_factor = discount_factor ;
parameters.interest_rate   = interest_rate   ;
parameters.initial_wealth  = initial_wealth  ;


% Creating a grid over assets
n_assets          = 200                     ;
assets            = linspace(0.01, 10, 200) ;
parameters.assets = assets                  ;

% Solving one period problem at t=2.  
continuation_values = zeros(1,n_assets) ;    % Setting continuation values V_3 = 0.
value_functions     = zeros(1,n_assets) ;    % First dimension is age, second is assets
for ia = 1:n_assets
    
    initial_guess = assets(ia);
    [ optimal_saving, optimal_utility ] = fminsearch(@(saving) solve_value(saving, assets(ia), continuation_values,...
                                                               assets, parameters), initial_guess); %#ok<ASGLU>
    value_functions(ia)                 = -1*optimal_utility ;
        
end
    
    
% Solving two period problem at t=1
continuation_values = value_functions      ;    % t=2 value function becomes t=1 continuation value.
% We will over-write the two-period value functions with the one-period
% value functions.  However, we will save the optimal savings this time.
savings_choices     = zeros( 1, n_assets ) ;
for ia = 1:n_assets
    
        initial_guess = assets(ia);
    [ optimal_saving, optimal_utility ] = fminsearch(@(saving) solve_value(saving, assets(ia), continuation_values,...
                                                               assets, parameters), initial_guess);
    value_functions(ia)                 = -1*optimal_utility ;
    savings_choices(ia)                 = optimal_saving     ;
    
end

figure
plot( assets, savings_choices)
xlabel('Assets')
ylabel('Optimal Saving')

figure
plot( assets, value_functions)
xlabel('Assets')
ylabel('Value Function')
        

end


function value = solve_value(saving, asset, continuation_values, assets, parameters)
    
    % Enforcing no-Ponzi condition
    if saving<0
        value = Inf;
        return
    end
    
    % Calculating consumption
    consumption = parameters.interest_rate*asset - saving;
    
    % Enforcing positive consumption
    if consumption<0
        value = Inf;
        return
    end
    
    % Calculating value
    value = ( 1/(1-parameters.crra) ) * ( consumption^( 1-parameters.crra ) )...
            + parameters.discount_factor * interp1(assets, continuation_values, saving, 'linear','extrap');
        
    value = -1*value;
end