% Loesung der Laplace Gleichung auf (0,1)^2 mit FDM
% uniforme Anordnung der Knoten 
% lexikographische Nummerierung der Knoten

clear all, close all, clc
numSteps = 6;
%*** definiere Volumenkraft f und Randbedingungen
global a;
a = 7;
u = @(x,y) sin(a*x).*sin(2*pi*y);
f = @(x,y) (a^2+4*pi^2)*sin(a*x).*sin(2*pi*y);
g = @(x,y,n) a*cos(a*x).*sin(2*pi*y)*n(1)...
        + 2*pi*cos(2*pi*y).*sin(a*x)*n(2);

% Definiere Gitterweite
h=1/3;

for k=1:numSteps
  % Definiere Anzahl der Punkte
  M=1/h-1;

  % Initialisiere Blockmatrix T (sparse)
  T = 1/h^2*spdiags([ -ones(M,1), 4*ones(M,1), -ones(M,1) ],[-1,0,1],M,M+1);

  % Assembliere Steifigkeitsmatrix A fuer Neumann-Bedingungen
  A = sparse(M^2+M,M^2+M);
  
  %*** TO  DO ***
  % setze Diagonalbl?cke
  for j =
    
  end
  % setze Nebendiagonalbl?cke
  for j = 
    
  end
  % setze Werte des Differenzensterns fuer Neumann-Knoten
  for j 
    
  end  
  %*

  % berechne rechte Seite fuer Neumann Randbedingungen
  s = linspace(0,1,M+2);
  [x,y] = meshgrid(s(2:end),s(2:end-1)); %*** Gitter freie Knoten
    
  %*** TO  DO ***
  % setze Vektor b
  b = 
  %*
  
  
  % Loesung des Gleichungssystems
  uh=A\b(:);

  % Graphische Darstellung
  U = zeros(M+2,M+2);
  U(2:M+1,2:M+2) = reshape(uh,M+1,M)';
  
  % Berechne Fehler auf dem groebsten Gitter
  idx = 2^(k-1)+1:2^(k-1):M+1;
  err(k) = max(max(abs( U(idx, idx)-u( x(idx-1,idx-1),y(idx-1,idx-1)) )));
  
  % Verfeinere Gitter
  h=h/2;

end

%*** TO  DO ***
figure(1)
% plotte Losung ueber dem Gitter

view(3);
title('Numerische Loesung der Laplace Gleichung')
xlabel('x')
ylabel('y')
zlabel('uh(x,y)');

figure(2)
% plotte Fehler (doppelt logarithmische Skala)

title('Maximaler absoluter Fehler auf dem grobsten Gitter')
xlabel('1/h')
ylabel('max|u-uh|')
loglogTriangle(1e1,1e2,1e-2,2,'l');
%*