function [T,Y] = heun(f,a,b,ya,m)
% [T,Y] = heun(f,a,b,ya,m)
% Heun's solution for y' = f(t,y) with y(a) = ya.
% f is the function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% ya is the initial condition, input.
% m is the number of steps, input.
% T is the vector of abscissas, output.
% Y is the vector of ordinates, output.
h = (b - a)/m;
T = zeros(1,m+1);
Y = zeros(1,m+1);
T(1) = a;
Y(1) = ya;
for j=1:m,
k1 = feval(f,T(j),Y(j));
p = Y(j) + h*k1;
T(j+1) = a + h*j;
k2 = feval(f,T(j+1),p);
Y(j+1) = Y(j) + h*(k1 + k2)/2;
end