function Lepuschitz

clear
clc

%----------- Einf. i. d. Numerik Übungsblatt 4 -----------
% Marlene Lepuschitz


% Anfangswerte:
N=25;
N1 = N+1;
N1N1 = N1*N1;
N4 = round(N/4);
V0 = [zeros(N4,1);ones(N1-2*N4,1);zeros(N4,1)];    
V0 = kron(V0,V0');
V0 = reshape(V0,N1N1,1);    
V=V0;
 
% erstellen von Matrix D

D0 = -eye(N1);
D1 = diag([4*ones(N,1);3])-diag(ones(N,1),+1)-diag(ones(N,1),-1);
D = kron(eye(N1),D1)+kron(diag(ones(N,1),+1),D0)+kron(diag(ones(N,1),-1),D0);
D2 = diag(ones(N1,1))-diag(ones(N,1),-1);

% erstellen von Matrix B
r = 1;
d = 1;
k = 1;
h = 1/N;
I = eye(N1N1);
F = kron(eye(N1),D2);
B = d/(h^2)*D+k/h*F+r*I;

% erstellen von Matrix A
M = 20;
tau= 1/M;
A = (I+tau*B);

%[x]=LoeseSystem(A,V);

% graphische Darstellung

V = V0;
for n=0:M
    V = LoeseSystem(A,V);
    surf(reshape(V,N1,N1)); axis([1 N1 1 N1 0 1]);
    view([1 1 1]); drawnow; pause(0.1);
end

end

function[x]= LoeseSystem(A,V)
   n   = size(A,1);
   d   = diag(A);
   x   = V;
   r   = x-A*x;
   tol = 1.0e-5;

   for k = 1:n
       z       = r./d;
       tau_new = z'*r;
       if (k==1)
           beta = 0;
           p    = z;
       else
           beta = tau_new/tau_old;
           p    = z+beta*p;
       end
       w     = A*p;
       a1    = tau_new / (p'*w);
       x_old = x;
       x     = x + a1 * p;
       if (norm(x-x_old) < tol*norm(x))
       return
       end
       r       = r-a1*w;
       tau_old = tau_new;


   end

end

