N  = 101;
N1 = N+1;
N4 = round(N/4);
h  = 1/N;
x  = linspace(0,1,N1);
M  = N;
t  = h;

r  = 1.0;
k  = 1.0;
d  = 1.0;

D = diag(2*ones(N1,1)) - diag(ones(N,1),1) - diag(ones(N,1),-1);
D(N1,N1) = 1;

F = diag(ones(N1,1)) - diag(ones(N,1),-1);

B = (d/h^2)*D + (k/h)*F + r*eye(N1);
A = eye(N1) + t*B;
[L,U,P] = lu(A);
if ((norm(P-eye(N1)) ~= 0) | ...
    (norm(L-tril(L)) ~= 0) | ...
    (norm(U-triu(U)) ~= 0))
   error('Permutation muss eingebaut werden...')
end

V = [zeros(N4,1);ones(N1-2*N4,1);zeros(N4,1)];
for n=1:M
    V = L\V;
    V = U\V;
    plot(x,V); axis([0 1 0 1]);
    drawnow; pause(0.1);
end
