function [] = Mein_AWP_Loeser(T,N)

close all; % Schließt alle Abbildungen

h = T/N;   % Schrittweite

fprintf('\nDie Schrittweite beträgt h = %1.10f\n',h);

% Initialisierungen
t = zeros(1,N+1);
y = ones(3,N+1);

% Parametervektor
% p = [0.1, 10, 100];             % -> Nicht-steifes System
p = [0.04, 10000, 30000000];      % -> Steifes System

% Knotenvektor
t = 0:h:T;

% Startwerte
y(:,1) = [1;0;0];


% Das klassische Runge-Kutta-Verfahren
% Stufe   s = 4
% Ordnung p = 4


for j=1:N
    
    k1 = f(t(j)     ,y(:,j)         ,p);
    k2 = f(t(j)+ h/2,y(:,j) + h/2*k1,p);
    k3 = f(t(j)+ h/2,y(:,j) + h/2*k2,p);
    k4 = f(t(j)+ h  ,y(:,j) + h  *k3,p);
    
    y(:,j+1) = y(:,j) + h*(1/6*(k1+k4) + 1/3*(k2+k3));
    
end


% Graphische Ausgabe
hold on;
plot(t,y(1,:),'b','LineWidth',2);
plot(t,y(2,:),'r','LineWidth',2);
plot(t,y(3,:),'k','LineWidth',2);
hold off;
box on;
legend('y1','y2','y3');

% Unterfunktion: Die rechte Seite der DGL
function [val] = f(t,y,p)

val = zeros(3,1);

val(1) = -p(1)*y(1) + p(2)*y(2)*y(3);
val(2) =  p(1)*y(1) - p(2)*y(2)*y(3) - p(3)*y(2)^2;
val(3) =                               p(3)*y(2)^2;
    



