% --------------------------------------------%
%      The Forward Euler Method for ODEs      %
%         Written by Prof. Dongwook Lee       %
%                AMS 20/20A                   %
% --------------------------------------------%


% ---------------------------------------------%
% Example 1:                                   %
% * dy/dt = f(t,y) = 3 - 2t - 0.5y             %
% * y(0) = 1                                   %
% ---------------------------------------------%

% [0] clear any previously defined values and
%     figures
clear all;
clf;

% [1] Set initial condition ====================
y0 = 1;



% [2] Set h ====================================
h = .25;



% [3] Set f(t,y) ===============================
f = inline('3 - 2*t - 0.5*y','t','y');



% [3] Temporal Evolution of ODE ================
% [3a] Set your clock to be zero initially
t = 0;

% [3b] Set your maximum timestep
n_step = 100;

% [3c] Set your maximum time based on the 
%      maximum timestep
t_max = 10;
n_max = 1000000;

% [3b] Temporal evolution of ODE until t=t_max
% Set your current y value y_n to be 
% the initial value to begin with
y = y0;

% Provide 1D arrays for solutions of length n_step
y_soln = zeros(n_step,1);
t_soln = zeros(n_step,1);


n = 1;
y_soln(n) = y0;

while (t<=t_max) & (n <= n_max);
    %for n=1:n_max;

        % Forward difference approximation
        y_new = y + h * f(t,y);
        
        % Update your clock
        t = t + h;
        n = n + 1;
        
        % Update your current value of y
        y = y_new;
        
        % Store your time and new solutions
        % in arrays for plot
        y_soln(n) = y;
        t_soln(n) = t;
    %end
    
end



% [4] Plot your numerical solution =============    
figure(1)
plot(t_soln(1:n), y_soln(1:n), 'r*')

% [5] Exact solution
y_exact = inline('14 - 4*t - 13*exp(-0.5*t)','t');

% [6] Overplot exact solution
hold on;
plot(t_soln(1:n), y_exact(t_soln(1:n)),'b');
grid on;
