function [x,varargout]=rkf45_g(f,x0,t0,t1,tol,xmax)
global previous_h

% Solving x'(t)=f(x(t)) 
% method: Runge-Kutta-Fehlberg 4th and 5th order

A=[        0,          0,          0,           0,         0;
         1/4,          0,          0,           0,         0;
        3/32,       9/32,          0,           0,         0;
   1932/2197, -7200/2197,  7296/2197,           0,         0;
     439/216,         -8,   3680/513,   -845/4104,         0;
       -8/27,          2, -3544/2565,   1859/4104,    -11/40];
b=[   16/135,          0, 6656/12825, 28561/56430,     -9/50,  2/55];
br=[   1/360,          0,  -128/4275, -2197/75240,      1/50,  2/55];
c=[ 0, 1/4, 3/8, 12/13, 1, 1/2];

n=length(x0);
nout=nargout-1; savepoints=0;
for j=1:nout , varargout(j)={NaN}; end
if nout > 0 , X=x0; T=t0; savepoints=1; end;
if nargin < 5 , tol=10^(-6); end;
if nargin < 6 , xmax=10^4; end;
k=zeros(n,6); h=previous_h; if isempty(h), h=0.0005;  end
t=t0; ready=0; x=x0;
hmin=0.00001; hmax=10; h=min([hmax,max(hmin,h),t1-t0]); 
if h<eps, return; end;
while not(ready) ,
   xo=x; 
   for j=1:6 , xs=xo+k(:,1:j-1)*A(j,1:j-1)'; k(:,j)=h*feval(f,t+c(j)*h,xs); end;
   err=norm(k*br')/h;  
   if err<tol , x=x+k*b'; t=t+h; if savepoints , X=[X,x]; T=[T,t]; end;end;
   h=min(hmax,max(min(4,0.85*(tol/err)^(1/4)),0.2)*h);
   if t1-t<h, previous_h=h; h=t1-t; 
      else , if h<hmin , 
      disp(sprintf('rkf45 failed, h=%0.5g t=%0.5g',h,t)), 
      f,x0,t0,t1,tol,x,err, x=NaN*x; return, end;
      if norm(x)>xmax , 
      disp(sprintf('rkf45 out of limits %0.5g',norm(x))),x=NaN*x; return, end;end;
   ready=(t>t1-tol); end;
if nout>0 , varargout(1)={X}; end;
if nout>1 , varargout(2)={T}; end;
