function [t,x]=heun(f,x0,t0,tf,h)
   t=[t0:h:tf];
   x=zeros(length(x0),length(t));
   x(:,1)=x0;
   for k=1:length(t)-1
       x(:,k+1)=heun_step(f,x(:,k),t(k),h);
   end 
endfunction

function xk1=heun_step(f,xk,tk,h)
    k1=f(xk,tk);
    k2=f(xk+h*k1,tk+h);
    xk1=xk+h/2*(k1+k2);
endfunction
