function [t,y] = euler(FunFcn,tspan,y0,n)
%
% This function produces an approximation to the solution of
% y' = f(t,y),  y(t0) = y0  at the times t = t0 ... tfinal
% at intervals of size h = (tfinal-t0)/n
% tspan = [t0, tfinal] is a vector of the initial and final times
% and y0 is the initial condition (it must be a column vector)
% h is the constant step size
% the output y is a matrix whose ith row contains an approximation
% to the solution [y_1(t0 + i*h), y_2(t0 + i*h), ..., y_m(t0+i*h)]
% i=1 ... n, where m is the number of components in the vector y
% the ouput t is a column vector containing the times t0 + h ... tfinal
% at which the approximate solution is computed
%
t0=tspan(1);
tfinal=tspan(2);
h= (tfinal-t0)/n;
n1 = n-1;
ynew = (y0 + h*feval(FunFcn, t0, y0));
y(1,:) = ynew';
tr(1) = t0 + h;
for i=1:n1
yold=ynew;
tr(i+1) = t0 + (i+1)*h;
ynew = yold + h*feval(FunFcn, tr(i), yold);
y(i+1,:) = ynew';
end
t = tr';