function UEB4_Bsp6_Glowatschnig
 
 
%---------------------
% Einführung in die Numerische Mathematik 1
% Alexander Glowatschnig
% Blatt 4 Übungsbeispiel 6
% --------------------
 
 
clear all;
clc
 
% Anfangswerte, Definitionen und Notationen
N=50;
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,N1);
D1 = diag([4*ones(N,1);3])-diag(ones(N,1),+1)-diag(ones(N,1),-1);
D = kron(eye(N1,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,N1),D2);
 
 
% konstruiere Matrix B:
d = 1;
k = 1;
r = 1;
h = 1/N;
I = eye(N1N1,N1N1);
B = d/(h^2)*D+k/h*F+r*I;
 
% erstelle Matrix A :
M = 20;
tau = 1/M;
A = (I+tau*B);
 
% graphische Darstellung
 
if ((d < 0) | ...
    (k < 0) | ...
    (r < 0))
    error('d, k, r duerfen nicht negative sein!')
end
for n=0:M
    if (k == 0)
        V = LoeseSystem(A,V);
    else
        V = LoeseSystem(A'*A,A'*V);
    end
    surf(reshape(V,N1,N1)); axis([1 N1 1 N1 0 1]);
    view([1 1 1]); drawnow; pause(0.1);
 end
 
end 
%Implementation von CG Verfahren
function x = LoeseSystem(A,b)
    n = size(A,1);
    d = diag(A);
    x = b;
    r = b-A*x;
    tol = 1.0e-4;  
     
    for k = 1:n
        z = r./d;
        t_new = z'*r;
        if (k==1) 
            be = 0;
            p = z;
        else 
            be = t_new/t_old; 
            p = z+be*p;
        end
        w  = A*p;
        al = t_new/(p'*w);
        x_old = x;
        x = x + al*p;
        if (norm(x-x_old) < tol*norm(x))
            return
        end    
        r = r-al*w;
        t_old = t_new;
        
    end
end
