function [t,x]=brk4(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)^4/24+(h*J)^3/6-(h*J)^2/2+h*J-eye(J);
       xktilde=rk4_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)^4/24+(h*J)^3/6-(h*J)^2/2+h*J-eye(J);
           xktilde=rk4_step(f,xold,t(k+1),-h);
           xnew=xold-inv(H)*(x(:,k)-xktilde);
       end
       x(:,k+1)=xnew;
   end 
endfunction

function xk1=rk4_step(f,xk,tk,h)
       k1=f(xk,tk);
       k2=f(xk+h/2*k1,tk+h/2);
       k3=f(xk+h/2*k2,tk+h/2);
       k4=f(xk+h*k3,tk+h);
       xk1=xk+h*(k1/6+k2/3+k3/3+k4/6);       
endfunction
