% Program 8.5 Finite difference solver for 2D Poisson equation 
%   with Dirichlet boundary conditions on a rectangle 
% Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps 
% Output: matrix w holding solution values 
% Example usage: w=poisson(0,1,1,2,4,4) 
function w=p6(xl,xr,yb,yt,M,N) 

P = 36.75;    % Power Wattage
K = 3.85;     % Thermal Conductivity (W/(cm*C))
H = 0.1;      % Convective Heat Transfer Coefficient (W/(cm^2*C))
L = 2;        % Power connection is 2 cm long
delta = 0.1;  % Width of fin is .1 cm

% Input Function is 0
f=@(x,y) 0;

m=M+1;           % Step size along x
n=N+1;           % Step size along y
mn=m*n;          % Total steps
A=zeros(mn,mn);  % Equation Matrix
b=zeros(mn,1);   % b vector for equation matrix

h=(xr-xl)/M;  % x Step Size
k=(yt-yb)/N;  % y Step Size
h2=h^2;       % Square of x step size
k2=k^2;       % Square of y step size

x=xl+(0:M)*h; % set mesh values 
y=yb+(0:N)*k; 

for i=2:M % interior points 
  for j=2:N 

    eqNumber = i + (j - 1)*m;  % Current Equation

    A(eqNumber, i - 1+(j-1)*m) = 1/h2;
    A(eqNumber, i + 1+(j-1)*m) = 1/h2; 
    A(eqNumber, i+(j-1)*m)= -2/h2 - 2/k2 - (2*H)/(K*delta);
    A(eqNumber, i+(j-2)*m) = 1/k2;
    A(eqNumber, i+j*m) = 1/k2; 
  end 
end 

for i=1:m % bottom and top boundary points 
  
  % Top
  j=n;
  A(i+(j-1)*m, i+(n-1)*m) = -3/(2*k) + H/K;
  A(i+(j-1)*m, i+(n-2)*m) = 2/k;
  A(i+(j-1)*m, i+(n-3)*m) = -1/(2*k);
  
  % Bottom
  j=1;
  A(i+(j-1)*m, i+(1-1)*m) = -3/(2*k) + H/K;
  A(i+(j-1)*m, i+(2-1)*m) = 2/k;
  A(i+(j-1)*m, i+(3-1)*m) = -1/(2*k);
end 

for j=2:n-1 % left and right boundary points 

  % Left Side
  i=1;
  if (j <= N/2.0)
     A(i+(j-1)*m, 1+(j-1)*m) = -3/(2*h); 
     A(i+(j-1)*m, 2+(j-1)*m) = 2/h; 
     A(i+(j-1)*m, 3+(j-1)*m) = -1/(2*h);
     b(i+(j-1)*m) = -P/(L*delta*K);  % Power
  else
     A(i+(j-1)*m, 1+(j-1)*m) = -3/(2*h) + H/K; 
     A(i+(j-1)*m, 2+(j-1)*m) = 2/h; 
     A(i+(j-1)*m, 3+(j-1)*m) = -1/(2*h);
  end
  
  % Right Side
  i=m;
  A(i+(j-1)*m, m+(j-1)*m) = -3/(2*h) + H/K;
  A(i+(j-1)*m, m-1+(j-1)*m) = 2/h;
  A(i+(j-1)*m, m-2+(j-1)*m) = -1/(2*h);
end 

v=A\b; % solve for solution in v labeling 
w=reshape(v(1:mn),m,n); %translate from v to w 
mesh(x,y,w') 

% Problem 6 Invokation:
% M = p6(0,4,0,4,10,10)
% M = p6(0,4,0,4,30,30)