function [y,t]=heun(f,a,b,y0,n)
h=(b-a)/n; 
t=a:h:b;
y=zeros(size(t)); 
y(1)=y0;
for k=1:n
yp= y(k)+h*f(t(k),y(k))       
y(k+1)=y(k)+(h/2)*(f(t(k),y(k))+f(t(k+1),yp));
end
t
y