function [t,x]=backheun(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)=backheun_step(f,x(:,k),t(k),h);
   end 
endfunction

function xk1=backheun_step(f,xk,tk,h)
   xk1=xk;//initial guess
   xk_l=heun_step(f,xk1,tk+h,-h);
   it=0;
   while norm(xk_l-xk)>1e-6*(norm(xk)+1e-3)&it<10
       J=jacobiant(f,xk1,tk+h);
       invH=pinv(eye(J)-J*h+(J*h)^2/2);
       xk1=xk1-invH*(xk_l-xk)
       xk_l=heun_step(f,xk1,tk+h,-h);
       it=it+1;
    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

function J=jacobiant(f,x,t)
 //estimates the Jacobian df/dx evaluated at (x,t) 
  n=length(x);
  J=zeros(n,n);
  fx=f(x,t);
  for i=1:n
    x1=x;
    dx=0.00001*x1(i)+1e-10;
    x1(i)=x1(i)+dx;
    J(:,i)=(f(x1,t)-fx)/dx;
  end
endfunction  