function [ x, y ] = euler ( f, xRange, yInitial, numSteps )
% [ x, y ] = euler ( f, xRange, yInitial, numSteps ) uses
% Euler's explicit method to solve a system of first-order ODEs
% y'=f(x,y).
% f = name of an m-file with signature
% yprime = f(x,y)
% to compute the right side of the ODE as a row vector
% xRange = [x1,x2] where the solution is sought on x1<=x<=x2
% yInitial = column vector of initial values for y at x1
% numSteps = number of equally-sized steps to take from x1 to x2
% x = row vector of values of x
% y = matrix whose n-th row is the approximate solution at x(n).
% M. M. Sussman
% $Id: euler.m,v 1.1 2007-12-25 18:45:45 mike Exp $
x=zeros(numSteps+1,1);
x(1) = xRange(1);
h = ( xRange(2) - xRange(1) ) / numSteps;
y(1,:) = transpose(yInitial);
for k = 1 : numSteps
x(k+1) = x(k) + h;
y(k+1,:) = y(k,:) + h * transpose(feval ( f, x(k), y(k,:) ));
end