function  [x, y] = eulers(f, ab, y0, N)
% Compute approximation of the solution of the initial value 
% problem  y' = f(x,y), y(a) = y0  on the interval ab = [a,b].
% Euler's method with step length  h = (b-a)/N .

% Version 11.12.2003.  INCBOX 

a = ab(1);  b = ab(2);  
h = (b - a)/N;
x = linspace(a,b,N+1);              % grid
y = zeros(size(x));  y(1) = y0;     % initialize  y
for  n = 1 : N
  y(n+1) = y(n) + h*feval(f, x(n),y(n));
end