function [u,A] = laplacefd(xx,top,bottom,yy,left,right)
%  Finite difference solution of Laplace's equation in a square
%  input parameters
%       xx   grid points for x coords  (row vector)
%      top   potential on top edge     (row vector)
%   bottom   potential on bottom edge  (row vector)
%       yy   grid points for y coords  (column vector)
%     left   potential on left edge    (column vector)
%    right   potential on right edge   (column vector)
%  output parameters
%        u   solution (matrix of all grid point values)
%        A   coefficient matrix
n=length(xx)-1;
nn=(n-1)*(n-1); h=xx(end)/n;
fprintf('subdivision parameter : %g\n',h)
% set up interior numbering
nint=reshape(1:nn,n-1,n-1);
zerorow=zeros(1,n+1); zerocol=zeros(n-1,1);
G=[zerorow; [zerocol,nint,zerocol]; zerorow];
% set up finite difference matrix using matlab function delsq
A=delsq(G);
% set up the RHS vector
b=zeros(nn,1);
b(1:n-1)=b(1:n-1)+left(2:n);
b(nn-n+2:nn)= b(nn-n+2:nn)+right(2:n);
b(1:n-1:nn)= b(1:n-1:nn)+top(2:n)';
b(n-1:n-1:nn)=b(n-1:n-1:nn)+bottom(2:n)';
bint=reshape(b,n-1,n-1);
% solve using matlab sparse solver
x=A\b;
xint=reshape(x,n-1,n-1);
% add boundary values and plot the solution
u=[top; [left(2:n),xint,right(2:n)]; bottom];
fdsolplot(xx,u)
return

function fdsolplot(xx,u)
[X,Y]=meshgrid(xx,xx);
subplot(121)
mesh(X,Y,flipud(u)), axis('square')
title('potential solution and grid')
subplot(122)
contour(X,Y, flipud(u),20), axis('square')
title('potential solution contours')
return