% Set RHS function, draw slope field & approx. solution with explicit Euler
% (EE) and implicit Euler (IE)
lambda = -1.;
f = @(t, y) lambda.*y;
tspan = [0., 8.5];
yspan = [-2., +2.];
clf;
draw_slope_field(f, tspan, yspan, 50);
hold on;
y0 = 1.;
[t_EE, y_EE] = euler_explicit(f, tspan, y0, 4);
plot(t_EE, y_EE, 'r-o');
[t_EE, y_EE] = euler_explicit(f, tspan, y0, 6);
plot(t_EE, y_EE, 'g-o');
[t_IE, y_IE] = euler_implicit(f, tspan, y0, 4);
plot(t_IE, y_IE, 'm-+');
legend('Slope field', ...
       'Explicit Euler (too large step size!)', ...
       'Explicit Euler (OK step size)', ...
       'Implicit Euler (No step size restriction!)');
hold off;


function hfig = draw_slope_field(f, tspan, yspan, Nsamp)
% draw_slope_field Calculate and draw the slope field of the right-hand side
% function of an ODE, i.e., dy/dt = f(t, y).
%
%   INPUT:  f     Function handle of t and y, calculating the right-hand side
%                 function of the dy/dt
%           tspan Time span over which to compute the vector field
%           yspan Solution span over which to compute the vector field
%           Nsamp Sampling number, i.e., how many "arrows" to draw in t and y
%
%   OUTPUT: hfig  Handle to the figure that is created

  t = linspace(tspan(1), tspan(2), Nsamp);
  y = linspace(yspan(1), yspan(2), Nsamp);
  [T, Y] = meshgrid(t, y);
  DT = ones(size(T));
  DYDT = f(T, Y);
  L = sqrt(DT.^2 + DYDT.^2);
  hfig = quiver(T, Y, DT./L, DYDT./L);
  xlabel('t');
  ylabel('y');
  axis tight;
end

function [t, y] = euler_explicit(f, tspan, y0, Nsteps)
  dt = (tspan(2) - tspan(1))/double(Nsteps);
  t = zeros(Nsteps + 1, 1);
  y = zeros(length(t), length(y0));
  t(1) = tspan(1);
  y(1) = y0;
  for j=1:Nsteps
    t(j+1) = t(j) + dt;
    y(j+1) = y(j) + dt*f(t(j), y(j));
  end
end

function [t, y] = euler_implicit(f, tspan, y0, Nsteps)
  dt = (tspan(2) - tspan(1))/double(Nsteps);
  t = zeros(Nsteps + 1, 1);
  y = zeros(length(t), length(y0));
  t(1) = tspan(1);
  y(1) = y0;
  for j=1:Nsteps
    t(j+1) = t(j) + dt;
    g = @(y_new) y_new - y(j) - dt*f(t(j+1), y_new);
    y(j+1) = fsolve(g, y(j), optimoptions('fsolve','Display','off'));
  end
end
