function [ t, y ] = euler ( dydt, tspan, y0, n )

%% euler approximates the solution to an ODE using Euler's method.
%
  tfirst = tspan(1);
  tlast = tspan(2);
  dt = ( tlast - tfirst ) / n;
  t = zeros ( n + 1, 1 );
  y = zeros ( n + 1, 1 );
  t(1) = tspan(1);
  y(1) = y0;

  for i = 1 : n
    t(i+1) = t(i) + dt;
    y(i+1) = y(i) + dt * dydt ( t(i), y(i) );
  end

  return
end

