clear all; clf % Löschen aller Variablen und Abbildungen aus dem Arbeitsspeicher

% Definition von Intergationsintervall und Variablen

g = 9.81;               % Erdbeschleunigung g in m/s^2
cd = 0.2 + 851/20000 ;  % Widerstandbeiwert plus letze drei Ziffern der
                        % Matrikelnummer geteilt durch 20000 in kg/m
m = 10;                 % Masse des Steines in kg

v_start = 0;

t0 = 0;                 % untere Grenze des Intervalls
tend = 10;              % obere Grenze des Intervalls
h = 1;
t = t0 : h : tend;
n = (tend - t0)/h;

%DGL als anonymous function
deltaV = @(V)   g-(cd/m)*V^2


% verbessertes Heun-Verfahren

V_hv = v_start;
Ea=1                        %Ea=1, für die Schlaufe
ind=0
while Ea>10^-7
    ind=ind+1
    k3 =deltaV(V_hv(ind));
    k4=deltaV(V_hv(ind)+(k3*h));
    V_hv(ind+1)=V_hv(ind)+(k3+k4)/2*h;
    
Ea=abs(((V_hv(ind+1)-V_hv(ind))/V_hv(ind+1)))

    if ind==100                 
    break
    end
  
end


% d) exakte Lösung der Fallgeschwindigkeit

V_exakt = sqrt(g*m/cd)*tanh(t*sqrt(g*cd/m))

%figure
plot(t,V_hv)    % k-- heißt K= schwarz, -- ist gestrichelt
grid on                            

legend('Fallgeschwindigkeit nach Heun in m/s','Fallgeschwindigkeit exakt in m/s')