function adaptiv_np(tn,t0,y0,h,tol)

% Die Funktion adaptiv_np loest ein AWP basierend auf dem Verfahren von
% Kutta und der verbesserten Polygonzugmethode mit einer adaptiven 
% Schrittweitensteuerung 
%
% Input: tn   Endzeit
%        t0   Startzeit
%        y0   Anfangswerte
%        h    Startschrittweite
%        tol  Toleranz zu Schrittweitenbestimmung
%
% Beispiele fuer Funktionsaufrufe:
%
% adaptiv_np(6.1,0,[1.2;0;0;-1.049357509830350],0.1,0.001)
% -> liefert periodische Loesung
%
% adaptiv_np(17.1,0,[0.994;0;0;-2.0015851],0.1,0.001)
% -> liefert chaotische Loesung

% Initialisierung
y(:,1)=y0(:);
t=t0;
T=t;
H=[];

% Integrationsschleife
while t<tn
    flag = 1;
    while flag
        % Hier soll das Programm entsprechend geaendert werden
        k1 = f(t,y(:,end));
        k2 = f(t+h/2,y(:,end)+h/2*k1);
        k3 = f(t+h,y(:,end)+h*(-k1+2*k2));
        dk = norm(h*(k1-2*k2+k3)/6);
        % Schrittweitensteuerung
        if dk > 1.2 * tol
            h = 0.9 * h;
        elseif dk < 0.8 *tol
            h = 1.1 * h;
        else
            flag = 0;
        end
    end

    y = [y,y(:,end)+h*k2];
    t=t+h;
    T=[T,t];
    H=[H,h];
end

% Grafische Ausgabe
figure(1)
semilogy([T(1:end-1);T(2:end)],[H;H],'g-');
title('Schrittweite');
print -depsc2 adaptiv_vp_fig01
figure(2)
plot(y(1,:),y(3,:),'r-');
title('Phasendiagramm');
print -depsc2 adaptiv_vp_fig02

end

% problemabhaengige Funktion f
function wert = f(t,y)
mu = 1/82.45;
z1 = ((y(1)+mu)^2+y(3)^2)^(3/2);
z2 = ((y(1)-1+mu)^2+y(3)^2)^(3/2);
wert = [y(2); y(1)+2*y(4)-(1-mu)*(y(1)+mu)/z1-mu*(y(1)-1+mu)/z2;
y(4); y(3)-2*y(2)-(1-mu)*y(3)/z1-mu*y(3)/z2];
end