% EX 7.16
% Bandmatrismetoden/Finita differensmetoden för linjärt problem 


% y'' - (2 - 2x + x^2)y' + y = 0, y(1) = 0, y(2) = 1
% Bestäm y(x) numeriskt. Beräkna y'(1) och int_1^2 y(x).
% Gör en tillförlitlighetsbedömning av resultaten.

figure(1)
clf
hold on
xlabel('X', 'FontSize', 14)
xlabel('Y', 'FontSize', 14)
axis([1 2 0 1.2])


disp('    h      yprim       Eyprim    Eyprimkvot=4     inty        Einty    Eintykvot=4')
step = 1;
Eyprim = [];
Einty = [];
yprim = [];
inty = [];
for n = [9 19 39 79]
	
	if (n == 9)
		col = 'ro';
	elseif (n == 19)
		col = 'go';
	elseif (n == 39)
		col = 'bo';
	else
		col = 'co';
	end
	
	% välj n punkter x_(i) = i*h med tillhörande värden y_(i)  
	% y_(0) = y(1), y_(n+1) = y(2)
	h = 1/(n + 1);
	y0 = 0;
	yend = 1;
	
	% Centraldifferens y_(i)'  = (y_(i+1) -  y_(i-1))/2h
	%                  y_(i)'' = (y_(i-1) - 2y_(i) + y_(i+1))/h^2
	% Matrisform: Ay = b, där A är en bandmatris (tridiagonal)
	% yprim(1) = (4y_(1) - y_(2))/2h  % andragrads Taylorutveckling
	% int(y) = (sum(y_(i)) + 1/2)/h   % Trapetsregeln
	x = 1 + (1:n)'*h;
	g = (2 - 2.*x + x.^2);
	A = spdiags([1/h^2 + g([2:end,end])/(2*h), (1 - 2/h^2)*ones(n,1), ...
		           1/h^2 - g([1,1:(end-1)])/(2*h)], -1:1, n, n);
	b = zeros(n, 1);
	b(n) = g(n)/(2*h) - 1/h^2;

	y = A\b;

	yprim = [yprim; (-3*y0 + 4*y(1) - y(2))/(2*h)];
	inty = [inty; h*(y0/2 + sum(y) + yend/2)];
	
	if (step > 1)
		Eyprim = [Eyprim; yprim(end-1) - yprim(end)];
		Einty = [Einty; inty(end-1) - inty(end)];
	end
	
	if (step > 2)
		fprintf(' %0.4f  %0.4e  %0.4e      %2.2f      %0.4e  %0.4e      %2.2f\n', h, ...
			yprim(end), Eyprim(end), Eyprim(end-1)/Eyprim(end), ...
			inty(end), Einty(end), Einty(end-1)/Einty(end));
	elseif (step > 1)
		fprintf(' %0.4f  %0.4e  %0.4e                %0.4e  %0.4e\n', h, ...
			yprim(end), Eyprim(end), inty(end), Einty(end));
	else
		fprintf(' %0.4f  %0.4e                            %0.4e\n', h, ...
			yprim(end), inty(end));
	end
	plot(x, y, col, 'LineWidth', 2)
	pause
	
	step = step*2;
end
	

% Körning:
%    h      yprim       Eyprim    Eyprimkvot=4     inty        Einty    Eintykvot=4
% 0.1000  6.1756e-01                            4.2515e-01
% 0.0500  6.1742e-01  1.4334e-04                4.2447e-01  6.8475e-04
% 0.0250  6.1736e-01  5.9760e-05      2.40      4.2430e-01  1.7133e-04      4.00
% 0.0125  6.1734e-01  1.7474e-05      3.42      4.2425e-01  4.2841e-05      4.00
