% Einführung in die Numerische Mathematik
% Aufgabe 6
% Gülnur Adanc


clc;
clear all;
close all;


% Anfangswerte
N = 10;
N1= N+1;
N4 = round (N/4);
V0 = [zeros(N4,1); ones(N1-2*N4,1); zeros(N4,1)];
nmax = 10;   % Anzahl der Zeitwerte

h = 1/N;
k = 0:1:N;
x = k*h;

% Erzeugung B:

d = 1;
k = 1;
tau = 1;
r = 1;
M = spdiags([2*ones(N,1);1],0,(N1),(N1)) + spdiags(-1*ones(N1),-1,(N1),(N1)) + spdiags(-1*ones(N1),1,(N1),(N1));
O = spdiags([1;ones(N-1,1);1],0,(N1),(N1)) + spdiags(-1*ones(N1),-1,(N1),(N1));

I = eye(N1);
B = (d/(h^2)) * M + (k/h) * O + (r * I);

% Erzeugung A
A = (eye(N1) + tau*B);

% Berechnung von V und graphische Darstellung
V= V0;   % V initialisieren

for n = 0:1:nmax
    
% LU-Zerlegung
[m,n] = size(A);
% if (m ~= n); error('Matrix A ist nicht quadratisch'); end
for k=1:n-1
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n) = A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);
end
L = eye(n,n) + tril(A,-1);U = triu(A);

% Vorwärts
n = size(L, 1);
V = zeros(1, n);

V(1) = V(1);
for i = 2 : n
V(i) = V(i) - sum(V(1:i-1) .* L(i, 1:i-1));
end

Vvorwaerts = V;

% Rückwärts
n = size(U, 1);
V = zeros(1, n);

V(n) = V(n) / U(n, n);
for i = n-1 : -1 : 1
V(i) = (V(i) - sum(V(i+1:n) .* U(i, i+1:n))) / U(i,i);
end

Vrueckwaerts = V;


% graphische Darstellung
plot(x,Vvorwaerts, x, Vrueckwaerts);
    drawnow; pause(0.1);
end



