function [tn,yn] = rkbs(f,t0,tend,y0,tol, h_start)
%
% Solve the ODE y'=f(t,y) with Bogacki-Shampine's method. 
% Input: f       : Function handle to the rhs. of the ODE. 
%        t0, tend: The integration interval. 
%        y0      : Initial value   
%        h       : Stepsize
% Output: tn     :
%         yn     : The solution at the points tn.

t = t0;    % Current solution
y = y0;

tn = t0;   % Store the solution
yn = y0;

h = h_start;
k1 = feval(f,t,y); 

while t < tend
    if t+h > tend
        h = tend-t;
    end
    
    % One step
    k2 = feval(f,t+0.5*h,y+0.5*h*k1);
    k3 = feval(f,t+0.75*h,y+0.75*h*k2);
    
    y_tmp = y+h*(2*k1+3*k2+4*k3)/9;
    k4 = f(t+h,y_tmp);
    
    % Error estimate
    errest = norm(h*(-5/72*k1+1/12*k2+1/9*k3-1/8*k4));
    if errest < tol
        % The step is successful
        y = y_tmp;
        t = t+h;
        k1 = k4;                   % FSAL property.
    
        % Store the solution
        tn = [tn,t];
        yn = [yn,y];
    end
    
    h = 0.9*(tol/errest)^(1/3)*h;
end


    
    



 


    



