function Holler

%% Einführung in die numerische Mathematik
%% Blatt 3
%% Beispiel 6

%Gernot Holler

%%----------------------------------------------

clear;clc;clf

%% Gegebene Daten
N=101;
d=0.09
h=1/N;
k=0.05
r=0.07
tau=h 

%% Bedingungen für tau, h >0    keine Pivotsuche notwendig falls
%% (d/h^2+k/h+r) != 1/tau

%% für d=0 r=0, k>0  -> reine Konvektion nach rechts
%% für d=0, k=0, r=0 -> wird an der Stelle abgebaut
%% für r=0, k=0, d>0 -> reine Diffusion, aber durch linken Rand begrenzt




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


%%Erstellen der gewünschten Matrix

A=-eye(N);
K=zeros(N,1);
C=zeros(1,N1);
I=eye(N1);
D=[ C ; A  K];
X=D+D'+2*I;
X(N1,N1)=1;
Y=D+I;

B= (d/h^2)*X+(k/h)*Y+r*I;
U=(I+tau*B);

%% Zerlegung in Matrizen L und U
for p=1:N
    
    p1=p+1;   
    
    mp(p)=U(p1,p)/U(p,p);
        
        for j=p:min(p+2,N1)  
            
            U(p1,j)=U(p1,j)-U(p,j)*mp(p);   
            
        end       
end
    
V=diag(mp);
L=[C;V K]+I;


%%Berechnung und Darstellung

nmax=N1;


V=V0; %Anfangsbedingungen

x = linspace(0,1,N1); %Vektor zur Darstellung

for n=1:nmax
    
    V=vorwaerts_gh(L,V);
    V=rueckwaerts_gh(U,V);
    plot(x,V);
    axis manual
    drawnow;pause(0.1);
    
end

end

function [y]=vorwaerts_gh(L,V)

N1=length(V);

y(1)=V(1);
for p=2:N1
    s=L(p,p-1)*y(p-1);
    y(p)=V(p)-s;
end

y=y';

end

function [x]=rueckwaerts_gh(U,V)

N1=length(V);

x(N1)=V(N1)/U(N1,N1);

for p=N1-1:-1:1
    
    s=U(p,p+1)*x(p+1);
    
    x(p)=(V(p)-s)/U(p,p);
    
end

x=x';

end
