function [tout, yout] = rungekuttaaut(F, RK, ab, y, tol)
%
% [tout, yout] = rungekutta(F, RK,ab, y0, k)
%
% INPUT:
% INPUT:
% F     - Function containing the ode
%         Call: yout = F(x,y)
%         x       - scalar x
%         y       - Solution column-vector.
%         yout    - Returned derivative column-vector; yout(i) = dy(i)/dx
% RK    - function with user-supplied RK-diagram.
%         [a,b,b2,c,ord]=RK()
%         a             - Matrix A
%         b             - vector b for order p
%         b2            - vector b for order q
%         c             - vector c
%         ord           - order q
% ab    - Initial value of a=ab(1), b=ab(2).
% y0    - Initial value column-vector.
% tol   - Tolerance
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% yout  - Returned solution, one solution column-vector per tout-value.
%
% Autor: Dieter Kaiser
%        Numerische Mathematik
%        Friedrich-Schiller-Universität Jena
%        http://users.minet.uni-jena.de/~kaiserd/
%        dieter.kaiser@uni-jena.de
% Date:  9-September-2013
[gamma,beta,beta2,rho,ord]=feval(RK);
y=y(:);
m=length(beta);dd=m-length(beta2);
n=length(y);
zu=1000;
tdim=zu;
yout=zeros(tdim,n);
i=1;
a=ab(1);b=ab(2);t=a;sh=sign(b-a);
tol=tol^1.1;
h=min(abs((b-a)/2),(tol^(1/(ord+1))/(norm(F(a,y))+1)))*sh;%h=(b-a)/2;
yout(1,:)=y';te=t+h;
tout=zeros(tdim,1);tout(1)=a;
zh=0;
while (b-te)*sign(h) >0
    te=t+h;
    if te == t,zh=zh+1;
        if zh == 8
            warning('on','AA:rkhsize')
            warning('AA:rkhsize','tout(end) ~= b')
            break
        end
        h=16*eps(t)*sh;
        te=t+h;
    end
    if (b-te)*sign(h) < 0,te=b;h=te-t;end    
    kj=zeros(n,m);
    hga=gamma*h;
    for j=1:m   
%Berechnung der kj - komplette Zeile
        kj(:,j)=F(t+h*rho(j),y+kj*hga(j,:)'); 
    end
%neues y für beide Verfahren
    yn=y+h*kj*beta';
    yn2=y+h*kj(:,1:m-dd)*beta2';
%Fehlerschätzung
    est=norm((yn-yn2))/(abs(h));
    %est=norm((yn-yn2),inf)/(abs(h));
    if ~isfinite(est),est=1/eps;end
    alf=min(.9*(tol/est)^(1/(ord+1)),2^10);
    h=alf*h;
    if est/tol < 1
        y=yn;t=te;
        i=i+1;
        if i>tdim, 
            yout=[yout;zeros(zu,n)];
            tout=[tout;zeros(zu,1)];
            tdim=tdim+zu;
        end
        yout(i,:)=y';
        tout(i)=t;
    else
    % schritt nicht erfolgreich - wiederholen   
    te=t;
    end
end
yout=yout(1:i,:);
tout=tout(1:i);