clear; clc; clf;
%% EINGABEWERTE

N = 150;

r = 0.8;
    % hohe reaktion 
    % -> kurve geht schneller gegen 0
k = 1.2;
    % konvektion 
    % -> kurve wandert nach rechts
d = 0.7;
    % hoher wert für diffusion 
    % -> kurve wird schneller flach
    
h = 1/N;
tau  = h;
nmax = floor(1/tau);

N1 = N+1;
N4 = round(N/4);
v  = [zeros(N4,1); ones(N1-2*N4,1); zeros(N4,1)];

%% DEFINIERE B, A
B =  d/h^2*(diag([2*ones(N,1);1]) ...
           - diag(ones(N,1),1) ... 
           - diag(ones(N,1),-1) ...
           ) ...
    + k/h*(eye(N1) - diag(ones(N,1),-1)) ...
    + r*eye(N+1);

A = eye(N1) + tau*B;

%% LU Zerlegung von B (tridiagonal)
for k = 1:N
    A(k+1, k)   = A(k+1,k)/A(k,k);
    A(k+1, k+1) = A(k+1,k+1) - A(k+1,k)*A(k,k+1);
end

%%

x   = zeros(N1,1);
int = linspace(0,1,N1);

for n = 1:nmax
    
    % Vorwärtssubstitution Lx = V
    x(1) = v(1);
    for i = 2:N1
        x(i) = v(i) - A(i,i-1)*x(i-1);
    end
    
    % Rückwärtssubstitution mit Uv = x
    v(N1) = x(N1)/A(N1,N1);
    for i = N:-1:1
        v(i) = (x(i) - A(i,i+1)*v(i+1)) / A(i,i);
    end
    
    % PLOT
    area(int,v);
    axis([-0.1 1.1 -0.1 1.1])
    axis equal
    drawnow; pause(0.1);
    
end
