function [U,err]=potentialcube(J)
% POTENTIALCUBE  Approximately solve a potential equation problem
%   u_xx + u_yy + u_zz = 0
% with particular Dirichlet bdry cond on the cube (0,1) x (0,1) x (0,1).
% Uses centered finite differences with equal spacing.  Solves linear
% system by using built-in sparse matrix solve (i.e. backslash).  Compares
% to exact solution.
% Form:
%     [U,err] = potential(J);
%   where
%           U = approx soln as (J+1)x(J+1)x(J+1) array; includes bdry vals
%         err = maximum error
%           J = number of grid spaces (in x, y, z)
% Example:  [U,err]=potentialcube(10);
% ELB 5/12/05
% See also: POTENTIAL2

dx=1/J;  dy=dx;  dz=dx;
% create A, b; solve
N=(J-1)^3;   b=zeros(N,1);
A=sparse(1:N,1:N,repmat(-6,1,N),N,N);  % fill in diagonal
for j=1:J-1
    for k=1:J-1
        for l=1:J-1
            q=QQ(J,j,k,l);
            % if neighbor is bdry w nonzero bdry cond, modify b
            if j==1, b(q)=-sin(pi*k*dy)*sin(pi*l*dz); end
            % otherwise put entry in A
            if j>1,  A(q,QQ(J,j-1,k,l))=1;  end
            if j<J-1,  A(q,QQ(J,j+1,k,l))=1;  end
            if k>1,  A(q,QQ(J,j,k-1,l))=1;  end
            if k<J-1,  A(q,QQ(J,j,k+1,l))=1;  end
            if l>1,  A(q,QQ(J,j,k,l-1))=1;  end
            if l<J-1,  A(q,QQ(J,j,k,l+1))=1;  end
        end
    end
end
% comment out for speed if J>15:  spy(A)
w=A\b;    % compute soln w (in column vector form)

% put result into 3d array U
U=zeros(J+1,J+1,J+1);  Uex=U;
% interior values
for j=1:J-1,  for k=1:J-1,  for l=1:J-1
    U(j+1,k+1,l+1)=w(QQ(J,j,k,l));
end,  end,  end
% nonzero boundary value:  u(0,x,y)=sin(pi x) sin(pi y)
for k=0:J,  for l=0:J
    U(1,k+1,l+1)=sin(pi*k*dy)*sin(pi*l*dz);
end, end

% compute exact soln on grid and error
for j=0:J,  for k=0:J,  for l=0:J
    Uex(j+1,k+1,l+1)=exact(j*dx,k*dy,l*dz);
end, end, end
err=max(max(max(abs(Uex-U))));

function q=QQ(J,j,k,l)
% go from grid (j,k,l) for U to global index q into vector w
q=(l-1)*(J-1)^2+(k-1)*(J-1)+j;

function u=exact(x,y,z)
% exact solution; requires scalar x, y, z (thus inefficient...)
s2p=sqrt(2)*pi;
xfact=( exp(s2p*x)-exp(2*s2p)*exp(-s2p*x) ) / (1-exp(2*s2p));
u=xfact*sin(pi*y)*sin(pi*z);