% Den här koden använder trapetsmetoden för att lösa ekvationen
%
%   y'(t) = t^2 * y,          y(0) = 1
%
% på tidsintervallet [0, 1]. Ekvationen är separabel, så vi kan
% hitta en exakt lösning analytiskt och därmed räkna ut felet
% exakt. (Baserat på Sauer Computer Problem 6.2.1)
%
% Joar Bagge 2016-11-11

clear

h = 0.1; % steglängd
y0 = 1; % initialvärde
t0 = 0; tend = 1; % tidsintervall

% Högerled till ekvationen y' = f(t,y)
f = @(t,y) t^2 * y;

N = (tend-t0) / h; % antal steg
t = t0:h:tend; % vektor med alla tidpunkter
y = zeros(1,length(t)); % lösningsvektor (fylls i loopen nedan)
y(1) = y0; % sätt initialvärdet!

% Trapetsmetoden
% Obs: de tre raderna inuti for-loopen är det enda som skiljer
% trapetsmetoden från Euler framåt!
for n=1:N
  k1 = f(t(n), y(n));
  k2 = f(t(n+1), y(n)+h*k1);
  y(n+1) = y(n) + h*(k1+k2)/2;
end

% Vi vet att den exakta lösningen är y=exp(t^3/3) (separabel ODE)
yexact = exp(t.^3 / 3); % räkna ut i varje tidpunkt
err = abs(y - yexact); % felet i varje tidpunkt

% Skriv ut en tabell med tidpunkterna, de numeriska värdena på y och felet
t = t'; y = y'; err = err'; % måste vara kolonnvektorer
table(t, y, err)

% Plotta lösningen och felet. Här ser man att det globala trunkeringsfelet
% växer med tiden.
figure(1)
plot(t, y)
xlabel('t'), ylabel('y'), title('Lösning')
figure(2)
plot(t, err, 'r')
xlabel('t'), ylabel('Fel'), title('Fel')
