function []=rungekutta(f,x0,T,x,tol)
%
% berechnet die Naeherungsloesung x_Delta
% des AWPs x'=f(x), x(0)=x0, 0 <= t <= T, mit exakter Loesung x
% mit dem expliziten Euler-Verfahren dem Verfahren von Runge 
% und  Runge-Kutta-4 zur aequidistanter Schrittweite bis zur Toleranz tol. 
% Plottet Naeherungsloesungen und Diskretisierungsfehler ueber f-Auswertungen.
%
% Parameter: 
% rechte Seite: f
% Startwert: x0
% Endzeit: T
% exakte Loesung: x
% Toleranz: tol
%
% Beispiele: rungekutta('fTan',0,1,'xTan',1e-3)
%            rungekutta('fTan',0,1,'xTan',1e-4)
%            rungekutta('fSqrt',1.e-3,1,'xSqrt',1e-1)
%

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

tau=T;
j=0;
n=1;
nRK=-1;

aEuler=1; aRunge=1; aRK4=1;

while aEuler | aRunge | aRK4 

  tau=tau/2;
  n=2*n;
  j=j+1;

  x_Euler(1)=x0;
  t_Euler(1)=0;
  errEuler(j)=0;
  evalEuler(j)=0;

  x_Runge(1)=x0;
  t_Runge(1)=0;
  errRunge(j)=0;
  evalRunge(j)=0;

  x_RK4(1)=x0;
  t_RK4(1)=0;
  errRK4(j)=0;
  evalRK4(j)=0;



  for k=1:n

    if aEuler
      t_Euler(k+1)=k*tau;
      x_Euler(k+1)=x_Euler(k)+tau*feval(f,x_Euler(k));
      errEuler(j)=max(errEuler(j),abs(feval(x,k*tau,x0)-x_Euler(k+1)));
      evalEuler(j)=evalEuler(j)+1;
    end

    if aRunge
      t_Runge(k+1)=k*tau;
      x_Runge(k+1)=x_Runge(k)+tau*feval(f,x_Runge(k)+0.5*tau*feval(f,x_Runge(k)));
      errRunge(j)=max(errRunge(j),abs(feval(x,k*tau,x0)-x_Runge(k+1)));
      evalRunge(j)=evalRunge(j)+2;
    end;

    if aRK4
      t_RK4(k+1)=k*tau;
      k1=feval(f,x_RK4(k));
      k2=feval(f,x_RK4(k)+0.5*tau*k1);
      k3=feval(f,x_RK4(k)+0.5*tau*k2);
      k4=feval(f,x_RK4(k)+1.0*tau*k3);
      sk=k1/6 + k2/3 + k3/3 + k4/6;
      x_RK4(k+1)=x_RK4(k)+tau*sk;
      errRK4(j)=max(errRK4(j),abs(feval(x,k*tau,x0)-x_RK4(k+1)));
      evalRK4(j)=evalRK4(j)+4;
    end;

  end;

  if (errEuler(j) < tol)&aEuler
    aEuler=0;
    tau_Euler=T/n
    j
    N_Euler=evalEuler(j)
    figure(1)
    nEuler=n;
    if nEuler< 65 cEuler='-g+'; else cEuler='-g'; end;
    if aRunge cRunge='-w'; else if nRunge< 65 cRunge='-c*'; else cRunge='-c'; end; end;
    if aRK4 cRK4='-w'; else if nRK4< 65 cRK4='-ro'; else cRK4='-r'; end; end;

    plot(t,x_Exakt,'k-',t_Euler,x_Euler,cEuler,t_Runge,x_Runge,cRunge,t_RK4,x_RK4,cRK4);
    xlabel('t')
    ylabel('solutions')
    legend('exact solution','Euler','Runge','Runge-Kutta-4');
    pause(1.e-4);
    hold off

    figure(2)
    cEuler='-g+';
    if aRunge cRunge='-w'; else cRunge='-c*'; end;
    if aRK4 cRK4='-w'; else cRK4='-ro'; end;
    loglog(evalEuler,errEuler,cEuler,evalRunge,errRunge,cRunge,evalRK4,errRK4,cRK4); 
    legend('Euler','Runge','RK4');
    xlabel('f-Auswertungen')
    ylabel('error')
    pause(1.e-4);
    hold off
  end;
  if (errRunge(j) < tol)&aRunge
    aRunge=0; 
    tau_Runge=T/n
    j
    N_Runge=evalRunge(j)
    figure(1)
    nRunge=n;
    if aEuler cEuler='-w'; else if nEuler< 65 cEuler='-g+'; else cEuler='-g'; end; end;
    if nRunge< 65 cRunge='-c*'; else cRunge='-c'; end;
    if aRK4 cRK4='-w'; else if nRK4< 65 cRK4='-ro'; else cRK4='-r'; end; end;

    plot(t,x_Exakt,'k-',t_Euler,x_Euler,cEuler,t_Runge,x_Runge,cRunge,t_RK4,x_RK4,cRK4);
    xlabel('t')
    ylabel('solutions')
    legend('exact solution','Euler','Runge','Runge-Kutta-4');
    pause(1.e-4);
    hold off

    figure(2)  
    if aEuler cEuler='-w'; else cEuler='-g+'; end;
    cRunge='-c*';
    if aRK4 cRK4='-w'; else cRK4='-ro'; end;
    loglog(evalEuler,errEuler,cEuler,evalRunge,errRunge,cRunge,evalRK4,errRK4,cRK4);
    legend('Euler','Runge','RK4');
    xlabel('f-Auswertungen')
    ylabel('error')
    pause(1.e-4);
    hold off
  end;

  if (errRK4(j) < tol)&aRK4 
    aRK4=0;  
    figure(1)
    j
    tau_RK4=T/n
    N_RK4=evalRK4(j)
    nRK4=n;
    if aRunge cRunge='-w'; else if nRunge< 65 cRunge='-c+'; else cRunge='-c'; end; end;
    if aEuler cEuler='-w'; else if nEuler< 65 cEuler='-g+'; else cEuler='-g'; end; end; 
    if nRK4< 65 cRK4='-ro'; else cRK4='-r'; end;
    plot(t,x_Exakt,'k-',t_Euler,x_Euler,cEuler,t_Runge,x_Runge,cRunge,t_RK4,x_RK4,cRK4);
    xlabel('t')
    ylabel('solutions')
    legend('exact solution','Euler','Runge','Runge-Kutta-4');
    pause(1.e-4);
    hold off

    figure(2)
    if aRunge cRunge='-w'; else cRunge='-c*'; end;
    if aEuler cEuler='-w'; else cEuler='-g+'; end;
    cRK4='-ro';
    loglog(evalEuler,errEuler,cEuler,evalRunge,errRunge,cRunge,evalRK4,errRK4,cRK4);
    legend('Euler','Runge','RK4');
    xlabel('f-Auswertungen')
    ylabel('error')
    pause(1.e-4);
    hold off  
   end;

end;

clear
