% 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
for i=2:9
    for j=2:9
        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
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];
DBCv = [0 0 0 0 0 0 0 0 0 0 20 20 20 20 20 20 20 20 20 20];
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))