function [U,err]=potential2(J,g,exactflag);
% POTENTIAL2  Approximately solve a potential equation problem
%   u_xx + u_yy = 0
% with a Dirichlet boundary condition  u|_{bdry}=g(x,y) on the square
% (0,1) x (0,1).  Uses centered finite differences with equal spacing
% in x and y.  Solves linear system by using built-in sparse matrix solve
% (i.e. backslash).  Compares to exact solution if available.
% Form:
%     [U,err] = potential(J,g,exactflag);
%   where 
%           U = approx soln as (J+1)x(J+1) matrix; includes bdry values
%         err = maximum error if computable (see exactflag)
%           J = number of grid spaces (in both x and y)
%           g = *inline* function of x and y which is value of u along bdry
%   exactflag = true if u(x,y)=g(x,y) is exact soln; false otherwise
% Example from complex vars; g(x,y) is real part of F(z) = e^{-2(z+1)^2}:
%    g=inline('exp(2.0*(y^2-(x+1)^2))*cos(4.0*(x+1)*y)','x','y')
%    [U,err]=potential2(50,g,true);
% Example, cont.; convergence and speed study:
%    JJ=[10 20 30 40 50 70 80 100]
%    for k=1:8, tic, [U,err(k)]=potential2(JJ(k),g,true); time(k)=toc, end
%    figure, loglog(1./JJ,err,'o'), xlabel('\Delta x'), ylabel('error')
%    perr=polyfit(log(1./JJ),log(err),1)
%    figure, loglog(JJ,time,'o'), xlabel J, ylabel('execution time')
% ELB 4/24/05; update of  potential.m  to make indexing match discussion
%   in class
% See also: DIFFUSION2, MINIMAL

dx=1/J;  dy=dx;

% create A, b; solve
N=(J-1)^2;   b=zeros(N,1);
A=sparse(1:N,1:N,repmat(-4,1,N),N,N);  % fill in diagonal with -4
for r=1:J-1
    for s=1:J-1
        k=kk(J,r,s);
        if r==J-1, b(k)=b(k)-g(1,s*dy);  % if right neighbor is bdry
        else, A(k,kk(J,r+1,s))=1; end  % otherwise, add to A
        if r==1, b(k)=b(k)-g(0,s*dy); % if left ...
        else, A(k,kk(J,r-1,s))=1; end
        if s==J-1, b(k)=b(k)-g(r*dx,1); % if upper ...
        else, A(k,kk(J,r,s+1))=1; end
        if s==1, b(k)=b(k)-g(r*dx,0); % if lower ...
        else, A(k,kk(J,r,s-1))=1; end
    end
end
w=A\b;    % compute soln w (in column vector form)

% put result into grid, with boundary values
U=zeros(J+1,J+1);
for r=1:J-1, for s=1:J-1,  U(r+1,s+1)=w(kk(J,r,s));  end, end  % interior
% lower and upper, the left and right boundaries
for r=0:J, U(r+1,1)=g(r*dx,0); U(r+1,J+1)=g(r*dx,1); end  
for s=1:J-1, U(1,s+1)=g(0,s*dy); U(J+1,s+1)=g(1,s*dy); end

% plot approx. soln; see "help mesh" on axis orientation
xx=0:dx:1;  yy=xx;  
figure(1), clf, mesh(xx,yy,U')
xlabel x, ylabel y, hidden off, title('approximate solution')

if exactflag  % compute and show error if possible; show matrix structure
    Uex=zeros(size(U));
    for r=0:J, for s=0:J, Uex(r+1,s+1)=g(r*dx,s*dy); end, end
    err=max(max(abs(Uex-U)));
    figure(2),  clf,  subplot(1,2,1),  mesh(xx,yy,abs(Uex-U)')
    xlabel x,  ylabel y,  hidden off,  title('absolute value of error')
    subplot(1,2,2),  spy(A),  title('structure of A')
    pos=get(gcf,'Position');  set(gcf,'Position',[pos(1) pos(2) 900 400])
else  % just show structure if exact solution not available
    err=NaN;  % note error is not computable
    figure(2),  clf,  spy(A),  title('structure of A')
end

function k=kk(J,r,s)  % go from grid (r,s) (for U) to vector index k
k=(s-1)*(J-1)+r;