function [t,y] = rk2(FunFcn,tspan,y0,n)
%
% This function uses Heun's method to produce
% 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;
k1 = feval(FunFcn, t0, y0);
ynew = y0 + .5*h*k1 + .5*h*feval(FunFcn,t0+h,y0+h*k1);
y(1,:) = ynew';
tr(1) = t0 + h;
for i=1:n1
yold=ynew;
tr(i+1) = t0 + (i+1)*h;
k1 = feval(FunFcn, tr(i), yold);
ynew = yold + .5*h*k1 + .5*h*feval(FunFcn, tr(i+1), yold+ h*k1);
y(i+1,:) = ynew';
end
t = tr';