% 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 # (j-1)*10 + i
nn = @(i,j) (j-1)*10 + i;

% 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; 3 2; 4 2; 5 2; 6 2; 7 2; 8 2; 9 2; ...
       2 3; 3 3; 4 3; 5 3; 6 3; 7 3; 8 3; 9 3; ...
   %   2 4; 3 4; 4 4; 5 4; 6 4; 7 4; 8 4; 9 4; ...
   %   2 5; 3 5; 4 5; 5 5; 6 5; 7 5; 8 5; 9 5; ...
       2 4;           5 4; 6 4; 7 4; 8 4; 9 4; ...
       2 5;           5 5; 6 5; 7 5; 8 5; 9 5; ...
       2 6; 3 6; 4 6; 5 6; 6 6; 7 6; 8 6; 9 6; ...
       2 7; 3 7; 4 7; 5 7; 6 7; 7 7; 8 7; 9 7;...
       2 8; 3 8; 4 8; 5 8; 6 8; 7 8; 8 8; 9 8;
       2 9; 3 9; 4 9; 5 9; 6 9; 7 9; 8 9; 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; 2 1; 3 1; 4 1; 5 1; 6 1; 7 1; 8 1; 9 1; 10 1; ...
       3 4; 4 4; 3 5; 4 5]; 
DBCv = [0 0 0 0 0 0 0 0 0 0 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 0 Neumann BCs
% this code is only good for the specific Neumann BCs in this problem
% which involve the x variable
NBC = [1 2; 10 2; 1 3; 10 3; 1 4; 10 4; 1 5; 10 5; 1 6; 10 6; 1 7; 10 7; ...
       1 8; 10 8; 1 9; 10 9];
nNBC = length(NBC);
for ii = 1:nNBC
    i = NBC(ii,1);
    j = NBC(ii,2);
    n = nn(i,j);
    A(n,n) = 1;
    if (i == 1) A(n,nn(2,j)) = -1;
    else A(n, nn(9,j)) = -1;
    end
   
    RHS(n) = 0;
end

% fill in equation matrix for grid points with Robin BCs
% this code is only good for the specific Robin BCs in this problem
% which involve the y variable
a = .1;
Tamb = 0;

for i = 1:10
    j = 10;
    n = nn(i,j);
    % the equation is T(i,10)(1 +dya)-T(i,9) = dyaTamb
    A(n,n) = 1 + dx*a;
    A(n, nn(i,9)) = -1;

    RHS(n) = dx*a*Tamb;
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))