% Löser y'(t) = f(t, y), y(0) = a, med Eulers metod, samt ritar ut
% riktningsfält och lösningskurva

% Mata in problemet
F = input('Ange högerled f(t, y), till exempel t*y^2: ', 's');
T = input('Ange sluttid > 0: ');
a = input('Ange begynnelsevillkor: ');
N = input('Ange antal steg: ');
eval(['f = @(t, y)', F, ';']) % tolka det inmatade högerledet som en funktion

% Lös med Eulers metod
h = T/N;                % steglängd
t = 0:h:T;              % sätt upp en vektor med tidpunkterna
y = zeros(N + 1, 1);    % sätt upp en vektor för lösningsvärdena
y(1) = a;               % begynnelsevillkoret vid starttiden t = 0
for i = 1:N;
    y(i + 1) = y(i) + h*f(t(i), y(i)); % ta ett steg
end

% Visualisera
ymax = max(y);          % största lösningsvärde 
ymin = min(y);          % minsta lösningsvärde
yspan = ymax - ymin;    % maximal variation i lösning (skala)
[tt, yy] = meshgrid(linspace(0, T, 18), linspace(ymin, ymax, 15)); % rutnät för riktningsfält
k = zeros(size(tt));    % skapa matris för riktningskoefficienter
for i = 1:numel(tt)
    k(i) = f(tt(i), yy(i)); % beräkna riktningskoefficienter
end
e = ones(size(tt)); % likformigt steg i t-led för riktningar
sc = sqrt((e/T).^2 + (k/1.2/yspan).^2); % skalningsmatris för riktningsfältet
quiver([tt, tt], [yy, yy], [e./sc, -e./sc], [k./sc, -k./sc], ...
       'ShowArrowHead', 'off', 'Marker', '.', 'MarkerSize', 10) % rita upp riktningsfältet
axis([0, T, ymin - 0.1*yspan - eps, ymax + 0.1*yspan + eps])
xlabel('t')
ylabel('y')
title('Riktningsfält (blå punkter och tangentlinjer) och lösning (ritas upp i rött)')
shg         % visa figur
hold on
pause(1.0)  % kort paus
for i = 1:N
    plot(t(i:i+1), y(i:i+1), 'r', 'LineWidth', 2)   % rita lösning bit för bit
    pause(0.2)                                      % ännu kortare paus
end
hold off

