%%
% ENM 7.16
%
% Lös randvärdesproblemet y'' - p(x)*y' + y = 0,
%                         y(1) = 0, y(2) = 1
% med finita differensmetoden.
%
% Joar Bagge 2016-11-29

clear

p = @(x) 2 - 2*x + x.^2;

% Randvillkor
xL = 1; yL = 0; % vänster rand
xR = 2; yR = 1; % höger rand

% Dela upp intervallet i n delintervall
n = 32;
x = linspace(xL, xR, n+1);
h = (xR-xL)/n; % delintervallens längd

% Approximera derivatorna med centraldifferenser...
% detta leder till ett linjärt system A*y=b.

%%
% Sätt upp systemmatrisen A
A = zeros(n+1, n+1);
% Översta och understa raden:
A(1,1) = 1; A(end,end) = 1
%%
% Mittenblocket:
for i=2:n
  A(i,i-1) = 1/h^2 + p(x(i))/(2*h);
  A(i,i) = 1 - 2/h^2;
  A(i,i+1) = 1/h^2 - p(x(i))/(2*h);
end
% (Istället för en for-loop går det att använda funktionen
% 'diag' eller 'spdiags' tillsammans med smart indexerande av A.)
A

%%
% Sätt upp högerledet b
b = zeros(n+1, 1);
b(end) = 1

%%
% Lös!
y = A\b
plot(x, y, '*-')
