function [] = solve_infinite_horizon()


% Assigning parameters
params.discount_factor = 0.96 ;
params.interest_rate   = 1.04 ;
params.wage            = 1.00 ;



% Discretizing the state space
asset_lower_bound = 0.001                                                    ;
asset_upper_bound = 100                                                      ;
n_assets          = 100                                                      ;
assets            = linspace(asset_lower_bound, asset_upper_bound, n_assets) ; 
params.assets     = assets                                                   ;

% Preallocating arrays
value_function      = zeros(1, n_assets) ;    % Also initial guess of value function.
savings             = zeros(1, n_assets) ;
new_value_function  = zeros(1, n_assets) ;

% Solving for the value function
tolerance = 10e-5 ;
iteration = 0     ;

while true
    iteration = iteration + 1;
    fprintf('Solving iteration %d.\n', iteration )
    
    for ia = 1:n_assets
        continuation_values = value_function ;
        asset               = assets(ia)     ;
        initial_guess       = assets(ia)     ;
        [optimal_saving, optimal_value] = fminsearch(@(saving) solve_value(saving, asset, continuation_values, params), initial_guess);
        new_value_function(ia) = -1*optimal_value;
        savings           (ia) = optimal_saving  ;
    end
    
    distance = max( abs( value_function - new_value_function ) ) ;
    fprintf('Current error is %0.5d.\n', distance )
    
    if distance<tolerance
        break
    end
    
    % Updating value function guess
    value_function = new_value_function;
    
end

figure
plot( assets, value_function )

figure
plot( assets, savings )



end

function lifetime_utility = solve_value(saving, asset, continuation_values, params)

% Testing Constraint
if ( saving<params.assets(1) ) || ( saving>params.assets(end) )
    lifetime_utility = inf;
    return
end

% Solving consumption and testing for positivity
consumption = params.interest_rate*asset - saving + params.wage ;

if ( consumption < 0 )
    lifetime_utility = inf;
    return
end

% Calculating continuation value
continuation_value = interp1( params.assets, continuation_values, saving, 'linear') ;

lifetime_utility = log( consumption ) + params.discount_factor * continuation_value ;
lifetime_utility = -1*lifetime_utility;



end