% function [y,t,hs]=rkf(a,b,y0,f,tol,hmax,hmin)
%
% Runge-Kutta-Fehlberg algorithm for solving IVP
% 
%    y'    = f(t,y),  a <= t <= b
%    y(a)  = y0
%
% see Burden and Faires, Algorithm 5.3
%
% inputs:
%  a,b   time interval to compute solution
%  y0    initial condition at t=a
%  f     function handle, for example
%
%        f = @(t,y) t^2 + y/2;
% 
%  tol   maximum error allowed
%  hmax  maximum step size allowed
%  hmin  minimum step size allowed
%
% outputs:
%  y     vector containing the approximation to y
%  t     vector of time history
%  hs    vector of time step history
function [y,t,hs]=rkf(a,b,y0,f,tol,hmax,hmin)
h=hmax; hs(1)=h; t(1)=a; y(1)=y0; flag=1; i=1;
% TO BE COMPLETED
