% ENM 7.16
% Denna kod skrevs på övningen 2017-10-11 i SF1514/18/19

% Problemet vi vill lösa är att bestämma y(x) så att
%
%   y'' - g(x)*y' + y = 0     för 1 <= x <= 2
%   y(1) = 0,
%   y(2) = 1.
%
% Vi gör detta med hjälp av finita differensmetoden. På övningen
% visade vi att problemet kan skrivas om som en matrisekvation
%
%   A*y = b.

clear

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

N = 8; % antal inre punkter
h = (2-1)/(N+1); % steglängd
x = 1:h:2; % alla x-punkter
xin = x(2:end-1); % alla inre x-punkter

% Skapa högerledet b
b = zeros(N, 1);
b(end) = -(1/h^2 - g(xin(N))/(2*h));

% Skapa matrisen A
A = zeros(N, N);
for k=1:N % alla rader i A
    A(k,k) = -2/h^2 + 1; % diagonalelementet
    if k < N % alla rader utom den sista
        A(k,k+1) = 1/h^2 - g(xin(k))/(2*h); % höger om diagonalen
    end
    if k > 1 % alla rader utom den första
        A(k,k-1) = 1/h^2 + g(xin(k))/(2*h); % vänster om diagonalen
    end
end
A, b

yin = A\b; % lösningen i de inre punkterna
y = [0; yin; 1]; % lösningen i alla punkter (lägg på randvillkor)
plot(x, y, '*-'), hold on
