function [ti,wi]=rkf(f,a,b,y0,tol,hmin,hmax)

%Runge-Kutta-Fehlberg method: See section 5.5.
%Solves y' = f(t,y), y(a) = y0, on the interval a <= t <= b.
%The unknown y may be a scalar or vector. (See section 5.9.)
%The code follows the implementation details from the book.
%The output ti is a column vector of values of t.
%The output wi is a column vector or matrix, whose rows are estimated values of y.

ti(1) = a;
wi(1,:) = y0(:);
h = hmax;
i = 1;
while true
  w = wi(i,:);
  t = ti(i);
  k1= h*f(t,w);
  k2= h*f(t+ 1/ 4*h,w+k1/4);
  k3= h*f(t+ 3/ 8*h,w+   3/  32*k1+   9/  32*k2);
  k4= h*f(t+12/13*h,w+1932/2197*k1-7200/2197*k2+7296/2197*k3);
  k5= h*f(t+      h,w+ 439/ 216*k1-   8     *k2+3680/ 513*k3- 845/ 4104*k4);
  k6= h*f(t+ 1/ 2*h,w-   8/  27*k1+   2     *k2-3544/2565*k3+1859/ 4104*k4-11/40*k5);
  % Compute 1/h times the difference between the 4th and 5th order estimates.
  diff = 1/h*abs(1/360*k1 - 128/4275*k3 - 2197/75240*k4 + 1/50*k5 + 2/55*k6);
  % Use the 4th order estimate and march forward one time-step if error is acceptable.
  if diff <= tol
    i = i+1;
    t = t+h;
    w = w + 25/216*k1 + 1408/2565*k3 + 2197/4104*k4 - 1/5*k5;
    ti(i,1) = t;
    wi(i,:) = w;
  end
  % Compute the factor by which h should be adjusted.
  % Assume that the true error is twice the 
  q = (tol/(2*diff))^(1/4);
  % Clamp this factor in the range [0.1, 4] to prevent large changes in step size.
  q = max(0.1, min(q,4));
  % The next step should have size q*h, unless that would
  % go outside the interval [a,b] or exceed the requested hmax.
  h = min([q*h, b-t, hmax]);
  if t >= b
    break;
  elseif h<hmin
    error('Minimum h exceeded');
  end
end