%6.Übung zur Numerischen Mathematik
%-Sparse-Systeme
clear; close all; clc;

%Festlegen der Rastergröße
M=input('Wie hoch soll das Raster sein?\nM = ');
N=input('Wie breit soll das Raster sein?\nN = ');

%Position des Spitzenwertes bestimmen
ny=ceil(1+(N-2)*rand(1));   
my=ceil(1+(M-2)*rand(1));

%Erzeugen der sparse-Matrix
L=(M-2)*(N-2);
A=sparse(L,L);

%Füllen der Matrix
k=0;
for n=2:N-1
    for m=2:M-1
        k=k+1; %Einzelindex inkrementieren
        
        a1=-4*(m+n);
        a2=m-1+n;
        a3=m+1+n;
        a4=m+n-1;
        a5=m+n+1;
        
        A_row=zeros(1,L);
        A_row(k)=a1;
        
        if(m~=2)
            A_row(k-1)=a2;
        end
        
        if(m~=M-1)
            A_row(k+1)=a3;
        end
        
        if(n~=2)
            A_row(k-M+2)=a4;
        end
        
        if(n~=N-1)
            A_row(k+M-2)=a5;
        end
        
        A(k,:)=A_row;
    end
end

%Einzelindex des Spitzenwertes (innerhalb der Matrix) festlegen
p=(ny-2)*(M-2)+(my-1);

%Gleichungssystem umschreiben, Zeile p löschen
A(p,:)=[];
b=-A(:,p);
A(:,p)=[];
z=A\b;

%Spitzenwert in z einfügen
z=[z(1:p-1); 1; z(p:(L-1))];

%Reshape des Vektors z zu einer Matrix der Rasterpunkte
Z=flipud((reshape(z,(M-2),(N-2))));

%Randnullen anfügen
Z=[Z zeros(size(Z,1),1)];   %linker Rand
Z=[zeros(size(Z,1),1) Z];   %rechter Rand
Z=[Z; zeros(1,size(Z,2))];  %unterer Rand
Z=[zeros(1,size(Z,2)); Z];  %oberer Rand
   
%Plot erzeugen
[X,Y]=meshgrid(1:N,M:-1:1);
surf(X,Y,Z);
title(['Spitzenwert bei (m,n)=(' num2str(my) ',' num2str(ny) ')']);
xlabel('n');
ylabel('m');