%-----------------------------------------------------------------
% Script 'Discretization Methods', Example 3.8
% (c) 08/2004 by A. Schmidt, IAM, Uni Stuttgart
%-----------------------------------------------------------------
clear all

% Variables
q_left = 0.0;                       % bound cond: heat flux at left bound
q_right = 0.0;                      % bound cond: heat flux at right bound
T_upper = 100.0;                    % bound cond: temperature at upper bound
T_lower = 50.0;                     % bound cond: temperature at lower bound
nx = 5;                             % number of sample points in x-direction
ny = 3;                             % number of sample points in y-direction
h = 1.0;                            % spacing of sample points 
                                    % (not needed for zero flux, but for plot)
k = 1.0;                            % thermal conductivity
                                    % (not needed for zero flux)

%-----------------------------------------------------------------
% Calculation of the System Matrix A: Az=b1-b2 [Eqs. (3.122)+(3.63), (3.124)]
% numbering of the sample points as given in the script
for i=1:1:ny                                  % number of 'blocks' in matrix A
  A((i-1)*nx+1,(i-1)*nx+1) = -4.0;            % first row in the 'block'
  A((i-1)*nx+1,(i-1)*nx+2) = 2.0;
  if i<ny
    A((i-1)*nx+1,i*nx+1) = 1.0;
  end
  if i>1
    A((i-1)*nx+1,(i-2)*nx+1) = 1.0;
  end

  for j=2:1:nx-1
    A((i-1)*nx+j,(i-1)*nx+j-1) = 1.0;         % middle of the 'block'
    A((i-1)*nx+j,(i-1)*nx+j) = -4.0;
    A((i-1)*nx+j,(i-1)*nx+j+1) = 1.0;
    if i<ny
      A((i-1)*nx+j,i*nx+j) = 1.0;
    end
    if i>1
      A((i-1)*nx+j,(i-2)*nx+j) = 1.0;
    end
  end

  A(i*nx,i*nx-1) = 2.0;                       % last row of the 'block'
  A(i*nx,i*nx) = -4.0;
  if i<ny
    A(i*nx,(i+1)*nx) = 1.0;
  end
  if i>1
    A(i*nx,(i-1)*nx) = 1.0;
  end
end

% Calculation of the right-hand side b1: Az=b=b1-b2
for i=1:nx:nx*ny
  b1(i) = -2.0*h/k*q_left;
end
for i=nx:nx:nx*ny
  b1(i) = -2.0*h/k*q_right;
end

% Calculation of the right-hand side b2: Az=b=b1-b2
for i=1:1:nx
  b2(i) = -T_upper;
end
for i=(ny-1)*nx+1:1:nx*ny
  b2(i) = -T_lower;    
end

b = (b1+b2)';                       % make b a column vector

%----------------------------
% Calculation of the vector of unknowns z (temperature): Az=b
z = A\b;

% Split temperatures in 2D array T(i,j) and add boundaries
for i=1:1:ny
  for j=1:1:nx
    T(i+1,j) = z((ny-i)*nx+j);
  end
end
T(1,1:nx) = T_lower;
T(ny+2,1:nx) = T_upper;

%-----------------------------------------------------------------
% Output
% axes
for i=1:1:nx
  x(i) = (i-1)*h;
end
for i=0:1:ny+1
  y(i+1) = i*h;
end

% plot
figure(1)
mesh(x,y,T)
xlabel('x / m')
ylabel('y / m')
zlabel('Temperature / °C')