function [yy, tt] = ode_rk54 (f, int, y0, tol, Nmax)

    % Coefficients for Prince-Dormand 5,4 pair:
    b(1,1)=0;
    b(2,1)=1/5;
    b(3,1)=3/40;         b(3,2)=9/40;
    b(4,1)=44/45;        b(4,2)=-56/15;         b(4,3)=32/9;
    b(5,1)=19372/6561;   b(5,2)=-25360/2187;    b(5,3)=64448/6561;
    b(5,4)=-212/729;  
    b(6,1)=9017/3168;    b(6,2)=-355/33;        b(6,3)=46732/5247;
    b(6,4)=49/176;       b(6,5)=-5103/18656;
    b(7,1)=35/384;       b(7,2)=0;              b(7,3)=500/1113;
    b(7,4)=125/192;      b(7,5)=-2187/6784;     b(7,6)=11/84;
    % 4th order scheme
    c4(1)=5179/57600;    c4(2)=0;               c4(3)=7571/16695;
    c4(4)=393/640;       c4(5)=-92097/339200;   c4(6)=187/2100;
    c4(7)=1/40;
    % 5th order scheme
    c5(1)=35/384;        c5(2)=0;               c5(3)=500/1113;
    c5(4)=125/192;       c5(5)=-2187/6784;      c5(6)=11/84;
    c5(7)=0;
    % Compute set of a-coefficients
    for i=1:7
      a(i) = sum(b(i,:));
    end

    t0   = int(1);
    t1   = int(2);
    h    = (t1-t0)/Nmax;      % initial step size
    hmax = (t1-t0)/10;
    hmin = (t1-t0)/1e9;
    y    = y0(:);             % Make sure we initialze with a column vector
    yy   = y;
    t    = t0;
    tt   = t0;
    step = 0;

    while t(length(t)) < t1

      step = step+1;
      if step >= Nmax
	 disp ('Maximum Number of steps reached before completion.');
	 return;
      end

      if h < hmin
	 disp ('The step size has become too small.  Singularity?');
	 return;
      end
      
      if t(length(t)) + h > t1
	 h = t1 - t(length(t));
      end

      % Since this is not a production code, we don't try to be clever.
      % For Prince-Dormand 5,4 the last stage of each step is actually
      % identical to the first stage of the next step, so we could
      % actually save one function evaluation...

      k(:,1) = feval (f, t, y);
      for j=2:7
	k(:,j) = feval (f, t + h*a(j), y + h*k(:,1:(j-1))*b(j,1:(j-1))');
      end

      y4 = y + h*k*c4;
      y5 = y + h*k*c5;

      error = norm (y4-y5);
      
      if error < tol
	 % Accept Step
	 y  = y5;
	 yy = [yy y];
	 t  = t+h;
	 tt = [tt t];
	 h  = min ([hmax, 2*h, 0.85*h*(tol/(eps+error))^(1/6)]);
      else % Reject Step
	 h = h/2;
      end

    end % Main loop

end
