function [t,x]=AB4(f,x0,t0,tf,h)
   t=[t0:h:tf];
   x=zeros(length(x0),length(t));
   fk=x;
   x(:,1)=x0;
  
   for k=1:3
       [x(:,k+1),fk(:,k)]=rk4_step_ab(f,x(:,k),t(k),h);
   end 
   
   for k=4:length(t)-1
      fk(:,k)=f(x(:,k),t(k));
      x(:,k+1)=x(:,k)+h/24*(55*fk(:,k)-59*fk(:,k-1)+37*fk(:,k-2)-9*fk(:,k-3));
   end 
endfunction


function [xk1,fk]=rk4_step_ab(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);
       fk=k1;       
endfunction
