function Adanc

clear all;
clc

% Anfangswerte, Definitionen und Notationen
N  = 20;
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;

% erstelle 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);
F  = kron(eye(N1),D2);


% konstruiere Matrix B:
d = 1;
k = 1;
r = 1;
h = 1/N;
I = eye(N1N1);
B = d/(h^2)*D+k/h*F+r*I;

% erstelle Matrix A :
M = 30;
tau = 1/M;
A = (I+tau*B);

% graphische Darstellung

for n=0:M
    V = LoeseSystem1(A'*A,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 = LoeseSystem1(A,b)
    n = size(A,1);
    d = diag(A);
    x = b;
    r = b-A*x;
    tol = 0.001;

    for l = 1:n
        m = r./d;
        t_new = m'*r;
        if (l==1)
            sigma = 0;
            p  = m;
        else
            sigma = t_new/t_old;
            p  = m + sigma*p;
        end
        v  = A*p;
        s = t_new / (p'*v);
        x_old = x;
        x  = x + s * p;
        if (norm(x-x_old) < tol*norm(x))
        return
        end
        r = r-s*v;
        t_old = t_new;


    end
end
