% 2-D steady state heat transfer
% model k*(d2T/dx2 + d2T/dy2) = 0 the Laplace equation
% 10x10 grid with Dirichlet BC and Neumann BC

T = zeros(10,10);     % array of temperature values
A = zeros(100,100);   % array of equation coefficients
RHS = zeros(100,1);   % array of equation right hand side values

dx=1;       % grid spacing
dy=1;

% assign node numbers, grid point (i,j) corresponds to node # (i-1)*10 + j
nn = @(i,j) (i-1)*10 + j;

% write equation matrix for interior grid points
% (T(i+1,j) - 2T(i,j) + T(i-1,j))/dx*dx +_(T(i,j+1) - 2T(i,j) + T(i,j-1))/dy*dy = 0

% list the interior grid points
INT = [2 2; 2 3; 2 4; 2 5; 2 6; 2 7; 2 8; 2 9; ...
       3 2; 3 3; 3 4; 3 5; 3 6; 3 7; 3 8; 3 9; ...
       4 2; 4 3; 4 4; 4 5; 4 6; 4 7; 4 8; 4 9; ...
       5 2; 5 3; 5 4; 5 5; 5 6; 5 7; 5 8; 5 9; ...
       6 2; 6 3; 6 4; 6 5; 6 6; 6 7; 6 8; 6 9; ...
       7 2; 7 3; 7 4; 7 7; 7 8; 7 9; ...
       8 2; 8 3; 8 4; 8 7; 8 8; 8 9
       9 2; 9 3; 9 4; 9 5; 9 6; 9 7; 9 8; 9 9];
nINT = length(INT);

for ii = 1:nINT
    i = INT(ii,1);
    j = INT(ii,2);
    n = nn(i,j);
    A(n, n) = -4;
    A(n, nn(i+1,j)) = 1;
    A(n, nn(i-1,j)) = 1;
    A(n, nn(i,j+1)) = 1;
    A(n, nn(i,j-1)) = 1;
    RHS(n) = 0;
end



% fill in equation matrix for grid points with Dirichlet BC
% list of grid points with Dirichlet BC
% this code is generic, i.e. it's good for Dirichlet boundary conditions
DBC = [ 1 1; 1 2; 1 3; 1 4; 1 5; 1 6; 1 7; 1 8; 1 9; 1 10; ...
        10 1; 10 2; 10 3; 10 4; 10 5; 10 6; 10 7; 10 8; 10 9; 10 10; ...
        7 5; 7 6; 8 5; 8 6]; 
DBCv = [0 0 0 0 0   0 0 0 0 0  5 5 5 5 5   5 5 5 5 5   10 10 10 10  ];
nDBC = length(DBC);
for ii = 1:nDBC
    i = DBC(ii,1);
    j = DBC(ii,2);
    n = nn(i,j);
    A(n,n) = 1;
    RHS(n) = DBCv(ii);
end

% fill in equation matrix for grid points with Newmann BC
% list of grid points with Dirichlet BC
% this code is only good for the specific Neumann BCs in this problem
% which involve the x variable
NBC = [2 1; 2 10; 3 1; 3 10; 4 1; 4 10; 5 1; 5 10; 6 1; 6 10; 7 1; 7 10; ...
       8 1; 8 10; 9 1; 9 10];
nNBC = length(NBC);
for ii = 1:nNBC
    i = NBC(ii,1);
    j = NBC(ii,2);
    n = nn(i,j);
    A(n,n) = 1;
    if (j == 1) A(n,nn(i,2)) = -1;
    else A(n, nn(i,9)) = -1;
    end
    RHS(n) = 0;
end

% let MATLAB compute the solution
sol = inv(A)*RHS;

% assemble the array for printing
for i = 1:10
    for j=1:10
        T(i,j) = sol(nn(i,j));
    end
end

surf(T)

% xi=1:dx:10;
% [xx,yy]=meshgrid(xi,xi);
% mesh(xx,yy,T(xi,xi))