% EX 3.10
% Olinjärt ekvationssystem 
% 1e5*x - 2e4*y + 1e4*x^3 = 0
% 1e6*x + 1e7*y + 5e6*y^2 = 3e6
clear

% Startgissning genom att ta bort alla olinjära termer och lösa det linjära
% ekvationssystemet
% 1e5*x - 2e4y = 0     
% 1e6*x + 1e7*y = 3e6           
ustart = [1e5 -2e4; 1e6 1e7]\[0; 3e6];
						
% Ekvationssystem, jacobian
F = @(u) [1e5*u(1) - 2e4*u(2) + 1e4*u(1)^3; ...
	        1e6*u(1) + 1e7*u(2) + 5e6*u(2)^2 - 3e6];
J_F = @(u) [1e5 + 3e4*u(1)^2, -2e4; ...
	          1e6,               1e7 + 1e7*u(2)];

% Fyra korrekta decimaler <=> absolutfel < 5e-5
% Men kör längre för att visa konvergensordning
disp(' n          u                      h            h./old_h.^2 = c (kvadr konv)')
fprintf(' 0  [%7.4f,%7.4f]\n', ustart(1), ustart(2));
u = ustart;
h = [1; 1];
n = 0;
while max(h) > 5e-9
	n = n + 1;
	old_h = h;
	h = J_F(u)\F(u);
	u = u - h;

	if n == 1
		fprintf('%2.0f  [%7.4f,%7.4f]  [%11.4e,%11.4e] \n', n, u(1), u(2), ...
			h(1), h(2));
	else
		c = h./old_h.^2;
		fprintf('%2.0f  [%7.4f,%7.4f]  [%11.4e,%11.4e]  [%7.4f,%7.4f]\n', ...
			n, u(1), u(2), h(1), h(2), c(1), c(2));

	end
	
end
		
% Körning:
% n          u                      h            h./old_h.^2 = c (kvadr konv)
% 0  [ 0.0588, 0.2941]
% 1  [ 0.0522, 0.2612]  [ 6.5961e-03, 3.2913e-02] 
% 2  [ 0.0521, 0.2608]  [ 8.5208e-05, 4.2269e-04]  [ 1.9584, 0.3902]
% 3  [ 0.0521, 0.2608]  [ 1.4051e-08, 6.9743e-08]  [ 1.9352, 0.3903]
% 4  [ 0.0521, 0.2608]  [ 3.8876e-16, 1.9267e-15]  [ 1.9692, 0.3961]