% ENM 7.4
% Denna kod skrevs på övningen 2017-10-05 i SF1514/18/19

clear

% Obs: vi kallar funktionsvariabeln för y här, fast i uppgiften är det h

for y0 = [0, 3, 10]
    % Gör framåt Euler för att hitta en approximativ lösning till ODE:n
    dt = 0.1/2; % steglängd
    t0 = 0;
    t_end = 40;
    num_steps = (t_end - t0)/dt;
    t = t0:dt:t_end; % skapa tidsvektor
    y = zeros(size(t)); % förbered lagring för lösningen till ODE:n
    y(1) = y0; % begynnelsevärde

    % Själva metoden:
    for n=1:num_steps
        y(n+1) = y(n) + dt*f(t(n), y(n));
    end
    plot(t,y)
    hold on
end


% Funktion för högerledet i y' = f(t,y)
function dhdt = f(t,h)
    A = 0.8; C = 0.1; g = 9.81;
    if h >= 0 && h <= 20
        Qin = 1 - h/20;
    else
        Qin = 0;
    end
    Qut = C*sqrt(2*g*h);
    dhdt = (Qin - Qut)/A;
end
