function ZP = potdiffS(N,bc)
%POTDIFFS  Solve potential eqn problem on a square using finite diffs.
%   Uses Matlab's sparse data structure, and so Matlab solves A\b
%   very quickly.  Allows input of boundary condition for y=1 side.
%   In calling "potdiffS(N,bc);", bc is a row vector of length N-1.
% Try "potdiffS(100,ones(1,99));" to solve problem 1.0.1.
% Try "potdiffS(100,jagged(.01:.01:.99,10,4,1.2));" for something different.
% (Ed Bueler; 2/7/02)
%
% Use finite differences to set up (N-1)^2 by (N-1)^2 matrix 
%    problem.  Unknown but pretty good scaling of work in N.
% Sets up matrix A using "numgrid" and "delsq"--this particular 
%    problem has simple enough geometry to allow their use.

% produce grid pattern in G
G=numgrid('S',N+1); % display G if desired with "G" or "spy(G)"

% ask delsq to compute the finite difference approximation to 
%  the matrix for  u_xx + u_yy
A=-delsq(G); % A is of sparse type because delsq is designed that way
% display A with spy if desired: figure(2); spy(A); 

b = zeros((N-1)^2,1); % not sparse is o.k. since Z=A\b won't be sparse
% fill in boundary values
b(N-1:N-1:(N-1)^2)=-bc;  % take advantage of Matlab!

% Solve A v = b.  Since A is sparse this is much faster and
%     uses less memory.  In fact, "A\b" is a very fancy command!
Z = A\b;

% turn Z back into grid coordinates for display
ZP = zeros(N+1,N+1);
ZP(G>0)=Z(G(G>0)); % sneaky way of putting values into grid
ZP(N+1,:)=[0 bc 0];
ZP=ZP';

% plot
dx=1/N; xp=0:dx:1; yp=xp;
figure(1); mesh(xp, yp, ZP);
xlabel('x axis'); ylabel('y axis'); axis([0 1 0 1 0 1.2]);
title('POTDIFFS: Solution of potential eqn problem on unit square.');