function Hrassnigg

clc;clf;clear;
%% beliebeige Werte (Achtung noch auszubessern)
N=100;
d=0.01; % d klein verursacht starke "Schwingung"
r=0.01; % r klein verursacht eine starke schwingung
k=1;  % k klein, lässt die Lösung sehr langsam konvergieren
tau=0.01; 

% die Parameter r,k d müssen so gewählt sein, dass A nicht singulär wird. 

%% Beginn code
h=1/N;
B1=zeros(N+1);
B2=eye(N+1);

% Matrix B erstellen
for i=1:N
    B1(i+1,i)=-1;
    B1(i,i+1)=-1;
    B1(i,i)=2;
    
    B2(i+1,i)=-1;
end

B1(N+1,N+1)=1;
I=eye(N+1);
 
B=d/h^2*B1+k/h*B2+r*I;
A=(I+tau*B);

if  det(A)==0
    error('Matrix ist singulär')
end

%% LU-Zerlegung von A:
L=eye(N+1);
Linv=L;
U=zeros(N+1);

for n=1:N
    M=eye(N+1);
    Minv=M;
    
        for i=n+1:N+1  %sollte stimmen
        M(i,n)=-A(i,n)/A(n,n);
        end
    
   
        for ij=n+1:N+1
        Minv(ij,n)=-M(ij,n);
        end
 

        A=M*A;
        U=A;
        L=L*Minv;
        Linv=M*Linv;
       
end

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

for i=1:N1
    x(i)=i*tau;
end
%% Vorwärts und Rückwärtssubstitution

V = V0;
for n=1:N1
V = vorwaerts(L,V); % Vorwaerts Substitution mit L
V = rueckwaerts(U,V); % Rueckwaers Substitution mit U
plot(x,V);
axis equal;
drawnow; pause(0.1);
end

end

function[u_k1]=vorwaerts(L,u_k)

sk=size(u_k,1);
u_k1=zeros([sk,1]);

u_k1(1,1)=u_k(1)/L(1,1);

for i=2:sk
    u_k1(i,1)=(u_k(i,1)-L(i,i-1)*u_k1(i-1,1))/L(i,i);
end

end

function[u_k1]=rueckwaerts(U,u_k)

sk=size(u_k,1);% ausgebessert? 2 statt1 

u_k1=zeros([sk,1]);
u_k1(sk)=u_k(sk)/U(sk,sk);

for i=sk-1:-1:1
    u_k1(i,1)=(u_k(i,1)-U(i,i+1)*u_k1(i+1,1))/(U(i,i));
end

end
