% euler_forward Euler approximation for ODE initial value problem
%               Forward Euler method
% [t, y] = euler_forward(f, t0, t1, y0, numSteps) returns a column vector of time
% points and a column vector of function values, 
% f = f(t, y) is the RHS, 
% t0 is the initial time, 
% t1 is the final time, 
% y0 is the initial function value, 
% numSteps is the number of steps of the numerical method
% 
% Lei Zhang, 11-09-2012
function [t, y] = euler_forward(f, t0, t1, y0, numSteps)

if nargin == 0
    test_euler_forward();
    return
end

% calculate time step dt
dt = (t1 - t0)/numSteps;

% initialize t and u as column vectors
t = linspace(t0, t1, numSteps+1);
y = zeros(size(t)); y(1) = y0;

% calculate u(t_n)
for k = 1:numSteps
    y(k+1) = y(k) + dt*f(t(k), y(k));
end

end

function test_euler_forward

choice = 1;
switch choice
    case 1
        f = @(t, u)pi*cos(pi*t);
    case 2
        f = @(t, u)(t + u);
    case 3
        f = @(t, u)exp(t);
    case 4
        f = @(t, u)-20*u + 20*t + 21;
    case 5
        f = @(t, u)u.^2.*(t-t.^3);
    case 6
        f = @(t, u)t^2 - u;
end
t0 = -1/2;
t1 = 3;
u0 = 1;

n = 2^16;
[t_ex, u_ex] = euler_forward(f, t0, t1, u0, n);
figure
plot(t_ex, u_ex)
hold

p = 3:10;
e = zeros(size(p));
for k = 1:length(p)
    n(k) = 2^p(k); 
    [t, u] = euler_forward(f, t0, t1, u0, n(k));
    if k==3
        plot(t, u, 'rx-')
    end
    e(k) = abs(u(end) - u_ex(end));
end
figure
loglog(n, e); 
P = polyfit(log(e), log(n), 1) 
hold on
loglog(n, 1./n, 'Color','r');
xlabel('n'); ylabel('error');
text(n(end), 1/n(end), 'n^{-1}');
end