% EX 7.8

% b)
% Beräkna IQ efter 100 dagar, q(100), mha RK4 med olika h
% q^2*q'' = 100, q(0) = 100, q'(0) = 0

% subst u1 = q, u2 = q':
% u1' = u2       , u1(0) = 100
% u2' = 100/u1^2 , u2(0) = 0

% Definiera funktion
f = @(t, u) [u(2); 100./u(1).^2];
 
figure(1)
clf
hold on
xlabel('Antal dagar i gymmet', 'FontSize', 14)
ylabel('IQ', 'FontSize', 14)
axis([0 100 100 160])

% RK4 med olika steglängd h / antal steg n
q100 = [];
disp('  h       q100        E_h      E_2h/E_h=16')
for n = [4 8 16 32 64]

	if (n == 4)
		col = 'r';
	elseif (n == 8) 
		col = 'g';
	elseif (n == 16)
		col = 'b';
	elseif (n == 32)
		col = 'c';
	else
		col = 'm';
	end

	% Steglängd given n lika långa steg
	h = 100/n;
	
	% Ta steg h i t-led, beräkna y
	t = 0;
	u = [100; 0];
	q = u(1);
	while t < 100
	
		% Ett steg i derivatans riktning
		k1 = h*f(t, u);
		k2 = h*f(t + h/2, u + k1/2);
		k3 = h*f(t + h/2, u + k2/2);
		k4 = h*f(t + h, u + k3);
		u = u + (k1 + 2*k2 + 2*k3 + k4)/6;
		q = [q, u(1)];
		t = t + h;
	
		plot([t - h, t], q([end-1, end]), col, 'LineWidth', 2)
		pause(0.1)
	end

	q100 = [q100, q(end)];
	if numel(q100) == 1
		fprintf('%5.2f  %10.6f\n', h, q100(end))
	elseif numel(q100) == 2
		fprintf('%5.2f  %10.6f  %0.4e\n', h, q100(end), ...
			(q100(end)-q100(end-1)))
	else
		fprintf('%5.2f  %10.6f  %0.4e  %7.2f\n', h, q100(end), ...
			(q100(end)-q100(end-1)), ...
			(q100(end-1)-q100(end-2))/(q100(end)-q100(end-1)))
	end
	pause
	
end

% Körning
%  h       q100        E_h      E_2h/E_h=16
%25.00  143.769296
%12.50  143.771238  1.9419e-03
% 6.25  143.771408  1.6998e-04    11.42
% 3.12  143.771420  1.2030e-05    14.13
% 1.56  143.771421  7.9383e-07    15.15



% c)
% omskrivning q' = sqrt(C - 200/q), q(0) = 100, q'(0) = 0.1
% 0.1 = sqrt(C - 2) => C = 2.01
% Lös med ode45, standardfelmarginaler 

% Definiera funktion
f = @(t, q) sqrt(2.01 - 200/q);

q0 = 100;
[t, q] = ode45(f, [0 100], q0);
q100 = q(end)

% Körning
%q100 =
%  151.6789