% 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=poisson2(0,1,1,2,4,4,5)

function w=poisson2(xl,xr,yb,yt,M,N,Power)
% inputs
d = .1;         % fin width in mm, delta
H = .005;       % W/cm^2 C, convective heat transfer coefficient
K = 1.68;       % W/cm C, thermal conductivity
L = 2;          % fin length in cm
P = Power;      % power in watts

m = M+1;
n = N+1;
mn = m*n;
h = (xr-xl)/M;
h2 = h^2;
k = (yt-yb)/N;
k2 = k^2;
x = xl + (0:M)*h;     % set mesh values
y = yb + (0:N)*k;
A = zeros(mn,mn);
b = zeros(mn,1);
F = 2*H/(K*d);

% interior points
for i=2:m-1
  for j=2:n-1
    A(i+(j-1)*m,i-1+(j-1)*m) = 1/h2;
    A(i+(j-1)*m,i+1+(j-1)*m) = 1/h2;
    A(i+(j-1)*m,i+(j-1)*m) = -((2/h2) + (2/k2)) - F;
    A(i+(j-1)*m,i+(j-2)*m) = 1/k2;
    A(i+(j-1)*m,i+j*m) = 1/k2;
    b(i+(j-1)*m) = 0;
  end
end

% bottom boundary
for i=1:m
   A(i,i) = -3/(2*k) + (H/K);
   A(i,i+m) = 4/(2*k);
   A(i,i+2*m) = -1/(2*k);
   b(i) = 0;
end;

% top boundary
for i=1:m
   A(i+(n-1)*m,i+(n-1)*m) = -(3/(2*k)) + (H/K);
   A(i+(n-1)*m,i+(n-2)*m) = (4/(2*k));
   A(i+(n-1)*m,i+(n-3)*m) = -(1/(2*k));
   b(i+(n-1)*m) = 0;
end

% left boundary
for j=2:n-1
    A(1+(j-1)*m,1+(j-1)*m) = -3/(2*h);
    A(1+(j-1)*m,2+(j-1)*m) = 4/(2*h);
    A(1+(j-1)*m,3+(j-1)*m) = -1/(2*h);
    b(1+(j-1)*m) = -P/(L*d*K);
end;

% right boundary
for j=2:n-1
    A(m+(j-1)*m,m+(j-1)*m) = -3/(2*h) + (H/K);
    A(m+(j-1)*m,m-1+(j-1)*m) = 4/(2*h);
    A(m+(j-1)*m,m-2+(j-1)*m) = -1/(2*h);
    b(j*m) =  0;
end;

% Solving
v = A\b;
w = reshape(v(1:mn),m,n);    % translating from v to w
mesh(x,y,w')
Error using poisson (line 13)
Not enough input arguments.