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


% y'' = - y + x, y(0) = 0, y(pi/2) = 1
% Analytisk lösning: y(x) = x + (1 - pi/2)sin(x)
% Bestäm y(x) numeriskt

figure(1)
clf
hold on
x = 0:(pi/200):(pi/2);
plot(x, x + (1 - pi/2)*sin(x), 'k', 'LineWidth', 2)
xlabel('X', 'FontSize', 14)
xlabel('Y', 'FontSize', 14)
pause

disp('    h     E_h(pi/8)     E_h(pi/4)     E_h(3*pi/8) mean(E_2h)/mean(E_h)=4')
step = 1;
x_E = pi/8*(1:3);
E = [];
for n = [3 7 15]
	
	if (n == 3)
		col = 'r*';
	elseif (n == 7)
		col = 'gx';
	else
		col = 'bo';
	end
	
	% välj n punkter x_(i) = i*h med tillhörande värden y_(i)  
	% y_(0) = y(0), y_(n+1) = y_(pi/2)
	h = (pi/2 - 0)/(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
	% i = 1         : y0 - 2y_(1) + y_(2) = h^2(- y_(1) + x)
	%                 (h^2 - 2)y_(1) + y_(2) = h^2x - y0
	% i = 2,..,(n-1): y_(i-1) - 2y_(i) + y_(i+1) = h^2(- y_(i) + x)
	%                 y_(i-1) + (h^2 - 2)y_(i) + y_(i-1) = h^2x
	% i = n         : yend - 2y_(n) + y_(n+1) = h^2(- y_(n) + x)
	%                 y_(n-1) + (h^2 - 2)y_(n) = h^2x - yend
	% Matrisform: Ay = b, där A är en bandmatris (tridiagonal)
	x = (1:n)'*h;
	A = spdiags([ones(n,1) (h^2 - 2)*ones(n,1) ones(n,1)], -1:1, n, n);
	b = h^2*x;
	b(1) = b(1) - y0;
	b(n) = b(n) - yend;

	y = A\b;

	E = [E; y(step:step:n)' - x_E - (1 - pi/2)*sin(x_E)];
	if (step > 1)
		fprintf(' %0.4f %0.6e %0.6e %0.6e      %2.2f\n', h, E(end, 1), E(end, 2), ...
			E(end, 3), mean(E(end - 1, :))/mean(E(end, :)));
	else
		fprintf(' %0.4f %0.6e %0.6e %0.6e\n', h, E(end, 1), E(end, 2), E(end, 3));
	end
	plot(x, y, col, 'LineWidth', 2)
	pause
	
	step = step*2;
end
	
