function [t,x]=bheun(f,x0,t0,tf,h)
   t=[t0:h:tf];
   x=zeros(length(x0),length(t));
   x(:,1)=x0;
   for k=1:length(t)-1
       xold=x(:,k);
       J=jacobian(f,xold,t(k+1));
       H=-(h*J)^2/2+h*J-eye(J);
       xktilde=heun_step(f,xold,t(k+1),-h);
       xnew=xold-inv(H)*(x(:,k)-xktilde);
       while norm(xnew-xold)>1e-8
           xold=xnew;
           J=jacobian(f,xold,t(k+1));
           H=-(h*J)^2/2+h*J-eye(J);
           xktilde=heun_step(f,xold,t(k+1),-h);
           xnew=xold-inv(H)*(x(:,k)-xktilde);
       end
       x(:,k+1)=xnew;
   end 
endfunction

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