function []=rungekuttaAdapt(f,x0,T,x,lTol)
%
% berechnet die Naeherungsloesung x_Delta
% des AWPs x'=f(x), x(0)=x0, 0 <= t <= T, mit exakter Loesung x
% mit dem Runge-Kutta-Fehlberg-Verfahren mit adaptiver Schrittweitensteuerung.
% Startschrittweite ist T/1000
% Plottet Naeherungsloesung, Schrittweiten und Fehler
% und gibt Diskretisierungsfehler und Anzahl der f-Auswertungen aus.
%
% Parameter:
% rechte Seite: f
% Startwert: x0
% Endzeit: T
% exakte Loesung: x
% lokale Toleranz und Startschrittweite: lTol
%
% Beispiel: rungekuttaAdapt('fTan',0,1,'xTan',1e-3)
%

t_Exakt=linspace(0,T,100);
x_Exakt=feval(x,t_Exakt,x0*ones(1,100));

tau_max=T/10;
rho=0.9; 

k=1;
x_RK4(1)=x0;
tau(1)=min([T/1000,tau_max]);

t_RK4(1)=0;

eval=0;
errRK4(1)=0;

while t_RK4(k)<T

   if k==1 k1=feval(f,x_RK4(k)); eval=eval+1;
   else if passed k1=k5; end; end;
   k2=feval(f,x_RK4(k)+0.5*tau(k)*k1);
   k3=feval(f,x_RK4(k)+0.5*tau(k)*k2);
   k4=feval(f,x_RK4(k)+1.0*tau(k)*k3);
   k5=feval(f,x_RK4(k)+tau(k)*(k1/6 + k2/3 + k3/3 + k4/6));
   eval=eval+4;

   epsk=tau(k)*abs(k4-k5)/6;

   if epsk<1.e-14 epsk=1.e-14; end;
   tauNew = min([2*tau(k),tau_max,tau(k)*(rho*lTol/epsk)^4]);
   tauNew = max([tau(k)/2,tauNew]);
   if tauNew <1.e-16
     'stepsize too small'
     break;
   end;

   if epsk<lTol
     t_RK4(k+1)=t_RK4(k)+tau(k);
     s_RK4=k1/6 + k2/3 + k3/3 + k4/6;
     x_RK4(k+1)=x_RK4(k)+tau(k)*s_RK4;

     errRK4(k+1)=abs(feval(x,t_RK4(k+1),x0)-x_RK4(k+1));

     tau(k+1)=min(T-t_RK4(k+1),tauNew);
     passed=1;
     k=k+1;
   else
     tau(k)=tauNew;
     passed=0;
   end;

end;


disp(sprintf('\n%d Funktionsauswertungen\n', eval));
disp(sprintf('Maximaler Fehler: %g\n', max(errRK4)));

figure(1)
clf
if eval<128 cRK4='r-o'; else cRK4='r-'; end;
plot(t_Exakt,x_Exakt,'k-')
xlabel('t');
ylabel('solutions');
%legend('exact solution','Runge-Kutta-4');
legend('exact solution');
hold on


mintau = min(tau);
if mintau<1e-10
   mintau = 1e-10;
end

figure(2)
clf
axis([0 T mintau max(tau)])
set(gca, 'YScale', 'log')
xlabel('t');
ylabel('step size');
hold on

minerr = min(errRK4);

if minerr<1e-20
   minerr = 1e-20;
end

figure(3)
clf;
axis([0 T minerr max(errRK4)])
set(gca, 'YScale', 'log')
xlabel('t');
ylabel('error');
hold on


n = size(t_RK4, 2);

for i=1:n-1

figure(1)
plot([t_RK4(i) t_RK4(i+1)], [x_RK4(i) x_RK4(i+1)], cRK4, 'EraseMode', 'none')

figure(2)
plot([t_RK4(i) t_RK4(i+1)], [tau(i) tau(i+1)], 'r-o', 'EraseMode', 'none')

figure(3)
plot([t_RK4(i) t_RK4(i+1)], [errRK4(i) errRK4(i+1)], 'r-o', 'EraseMode', 'none')

pause(0.2);
end


clear
