function [t,x]=rk4(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)=rk4_step(f,x(:,k),t(k),h);
   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/6*(k1+2*k2+2*k3+k4);
endfunction
