% Euler Framåt
% du/dt = -u^2 + t^2
% u(1) = 1
clear

% Definiera funktion
f = @(t, u) -u.^2 + t.^2;

% Rita riktningsfält
[tkoord, ukoord] = meshgrid(1:0.05:2, 1:0.05:2);
pilt = ones(size(tkoord));
pilu = f(tkoord, ukoord);

figure(1)
clf
hold on
quiver(tkoord, ukoord, pilt, pilu, 'k')
xlabel('t')
ylabel('u')
axis equal
axis tight
axis([1 2 1 2])

% "sann" lösning med ode45(f, [t_start t_end], y_start)
[true_t, true_u] = ode45(f, [1 2], 1);
plot(true_t, true_u, 'k', 'LineWidth', 2)
pause

% Euler framåt med olika steglängd
E = 1;
disp('   h       E_h   E_2h/E_h=2')
for h = [0.5 0.25 0.125 0.0625]

	if (h == 0.0625)
		col = 'r';
	elseif (h == 0.125)
		col = 'g';
	elseif (h == 0.25)
		col = 'b';
	elseif (h == 0.5)
		col = 'c';
	end

	
	% De t för vilka vi beräknar u
	all_t = 1:h:2;

	% Ta steg h i t-led, beräkna u
	t = 1;
	u = 1;
	all_u(1) = u;
	step = 1;
	while t < 2
	
		% Ett steg i derivatans riktning
		u = u + h*f(t, u);

		step = step + 1;
		t = t + h;
		all_u(step) = u;
	
		plot(all_t([step-1 step]), all_u([step-1 step]), col, 'LineWidth', 2)
		pause(0.1)
	end

	E = [E true_u(end) - u];
	fprintf('%0.4f  %9.3e  %6.3f\n', h, E(end), E(end-1)/E(end))
	pause
	
end
disp('Linjärt avtagande trunkeringsfel - noggrannhetsordning p = 1')
