%% Einfache Explizite Einschrittverfahren 
% Für eine explizite Differentialgleichung 1. Ordnung
%
% $$y' = f(x,y)$$
%
% mit Anfangsbedingung ist hier der Ablauf eines Einschritt-Verfahrens
% programmiert.
%
% Die Funktion f(x,y) und die Verfahrensfunktion F(x,y,h) sind als anonymous functions implementiert. 
%
% Version 2020: anonymous functions               CWB SS 2014,16,19,20

%% Die Funktion f definiert die Differentialgleichung
f = @(x,y) 1/(y+1)-x/4;
%%
% Zum Vergleich: Sie könnten $f$ so auch in einer eigenen Funktionsdatei
% oder als Subfunktion am Ende des Skript definieren (hier auskommentiert)
%%
%   function ystrich = f(x,y)
%      ystrich = 1/(y+1)-x/4;
%   end


%% Die Verfahrensfunktion F
% Die Verfahrensfunktion F gibt die aktuelle Fortschreit-Richtung an. Je
% nach Rechenverfahren ist sie unterschiedlich zu programmieren. Hier 
% eine Implementierung des Heun-Vefahrens 2. Ordnung als anonymous function
F = @(x,y,h) (f(x,y) + f(x+h, y+h*f(x,y)))/2;  

% F = @(x,y,h) f(x,y);  % Einfaches explizites Euler-Verfahren

%%
% Die Implementierung als anonymous-function-Einzeiler wird für etwas
% komplexere Verfahrensfunktionen grausam unübersichtlich.
% Programmiertechnisch übersichtlicher ist eine Subfunktion oder
% eigene Funktionsdatei. 
% Warum machen wir das dann nicht? 
%%
%
% # Subfunktionen in Skript-Dateien funktionieren in der alten
% MATLAB-Versio R2013a auf meinem Präsentations-Laptop nicht
% # Subfunktionen in Skripts dürfen erst ganz am Ende der Datei
% stehen; übersichtlicher ist es, wenn die Angabe zu Beginn steht.
%
% Zum Vergleich: so könnte eine Subfunktion aussehen (hier auskommentiert)
%%
%   function richtung = F(x,y,h)
%      k1 = f(x,y);
%      k2 = f(x+h, y+h*k1);
%      richtung  = (k1 + k2)/2;    % hier ist das Heun-Verfahren implementiert
%   end


%% Vorbereitung zur numerischen Lösung: Setze Anfangwerte, Intervall, Schrittweite
% Anfangswerte
x = 0;
y = 2;
%%
% Endwert des x-Intervalls, Schrittweite
x_end = 3;
h = 1/2;
%%
% Anzahl der Rechenschritte
n = round((x_end - x)/h);
%%
% Es kann sein, dass sich $x_end$ nicht genau in $n$ Schritten erreichen
% laesst - darum kümmern wir uns hier nicht.

%%
% Initialisiere die Lösungs-Vektoren
xLsg = zeros(1,n);
xLsg(1) = x;
yLsg = zeros(1,n);
yLsg(1) = y;
 
%% Einschritt-Verfahren: Schleife

for i=1:n
    y = y + h*F(x,y,h);         % berechne neues x und y
    x = x + h;                  % (Achtung: zuerst y, weil dafuer noch
    % alter x-Wert eingesetzt wird!)
    xLsg(i) = x;       	% speichere Werte in Loesungsvektor
    yLsg(i) = y;
end

%% Ausgabe: Plot und Endwert
plot(xLsg, yLsg,'-o');
fprintf('Endwert: (x, y ) = ( %g , %g )\n', x, y);%%





