%---------------------------------------------------
function [t,y,Fevals,rejects] = rk4adaptive(f,y0,a,b,tol)
%input: function handle f auf rechte Seite
%       Startvektor y_0
%       Startzeitpunkt a, Endzeitpunkt b
%       Toleranz tol 
%
%output: t = Vektor mit gewaehlten Knoten t_i
%        y = Vektor mit zugehoerigen Werten
%        Fevals = Anzahl Funktionsauswertungen
%        rejects = Anzahl 'step rejects'
%
    h = tol^(1/5);
    hmin = tol;
    rho = 0.8;
    eta = 2;
    p = 4;
    Fevals = 0;
    rejects = 0;
%    NMAX = ceil((b-a)/tol);
    NMAX = 100000; 
    nu = size(y0,1);
    t = zeros(1,NMAX);
    y = zeros(nu,NMAX);
    y(:,1) = y0;
    t(1) = a;
    i = 1;
    h = min(h,b-a); 
    noch_nicht_schluss = true; 
    while ( noch_nicht_schluss ) 
        ytilde = rk4step(f,y(:,i),t(i),h);
        y12 = rk4step(f,y(:,i),t(i),h/2);
        yhat = rk4step(f,y12,t(i)+h/2,h/2);
        Fevals = Fevals + 12;    %11 ginge auch ...
        EST = norm(ytilde-yhat) / ( 1 - 2^(-p) );
        if ( ( EST/h) <= tol || h <= hmin )
            % accept
            %disp(['accept ',num2str(t(i)),' ',num2str(h)]);
            y(:,i+1) = yhat;
            t(i+1) = t(i) + h;
            if ( t(i+1)==b )
               noch_nicht_schluss = false; 
            end
            h = max(hmin,min(eta*h,rho*(tol/EST*h^(p+1))^(1/p)));
            if ( t(i+1)+h >= b ) 
                h = b-t(i+1);
            end
            i = i+1;
        else
            %disp(['reject',num2str(t(i)),num2str(h)]);
            h = h/2;
            
            rejects = rejects + 1;
            Fevals = Fevals - 4;   %beim "step reject" auch Wiederverwendung denkbar...
        end
    end
    
    y = y(:,1:i);
    t = t(1:i);
return
end
%-----------------------------------------------------------------    
    function y = rk4step(f,y0,t,h)
        k1 = feval(f,t,y0);
        k2 = feval(f,t+h/2,y0+h/2*k1);
        k3 = feval(f,t+h/2,y0+h/2*k2);
        k4 = feval(f,t+h,y0+h*k3);
        y = y0 + h*( k1+2*k2+2*k3+k4 ) / 6;
   return
end
%-----------------------------------------------------------------