function U=minimal(J,g,TOL,f,U0);
% MINIMAL  Solves the minimal surface equation
%   div( a(grad u)  grad u )  + f = 0
% where a(grad u)=(1+u_x^2+u_y^2)^{-1/2}.  Requires boundary condition 
% g=g(x,y), input as an inline function or function handle.  Accepts 
% forcing term f(x,y), input as an inline function or function handle, 
% which must be written to take vector arguments.  Default f(x,y)=0.  (See 
% examples below.)
% Does naive iteration; the default values to start with, U0, is to compute
% solution of  div(grad u) = 0  with given boundary condition.  U0 may also
% be given by user.
% Form:
%   U=minimal(J,g,TOL,f,U0)
% where
%     U = approximate solution ( (J+1)x(J+1) matrix)
%     J = number of grid spaces (in both x and y)
%     g = g(x,y) = inline function or function handle; u on bdry
%   TOL = iteration tolerance (for max diff between iterates)
%    [f = f(x,y) = inline function or function handle; source;  OPTIONAL]
%   [U0 = (J+1)x(J+1) matrix which is first guess for soln; OPTIONAL] 
% Example I, with comparison to solution to linear eqn div(grad u)=0:
%   >> J=50;  g=inline('4*x.*y','x','y') 
%   >> U=minimal(J,g,1e-4);
%   >> figure, Ulin=diffusion2(J,ones(J,J-1),ones(J-1,J),zeros(J-1),g);
% Example II: USEMIN
% Example III: USEMIN_SMART
% ELB 4/24/05
% See also: DIFFUSION2, USEMIN, USEMIN_SMART

maxiter=300;
dx=1/J;  dy=dx;  x=0:dx:1;  y=x;  % x,y have length J+1
[XX,YY]=ndgrid(x(2:J),y(2:J));  % interior pts; for eval of forcing

% use f if given
if (nargin==4)|(nargin==5),  fmatrix=feval(f,XX,YY);
elseif nargin==3,  fmatrix=zeros(J-1,J-1);
else, error('must have 3, 4, or 5 input arguments'), end  

% compute U, or use input U, to get started
if nargin==5
    U=U0; % set U if known
else
    aLEFT=ones(J,J-1); aDOWN=ones(J-1,J);  % default
    U=diffusion2(J,aLEFT,aDOWN,fmatrix,g);
end 

for m=1:maxiter
    dxLEFT=(U(2:J+1,2:J)-U(1:J,2:J))/dx;
    dyLEFT=(  ( U(1:J,3:J+1)+U(2:J+1,3:J+1) )/2 - ...
              ( U(1:J,1:J-1)+U(2:J+1,1:J-1) )/2  )/(2*dy);
    dxDOWN=(  ( U(3:J+1,2:J+1)+U(3:J+1,1:J) )/2 - ...
              ( U(1:J-1,2:J+1)+U(1:J-1,1:J) )/2  )/(2*dx);
    dyDOWN=(U(2:J,2:J+1)-U(2:J,1:J))/dy;
    aLEFT=( 1+dxLEFT.^2+dyLEFT.^2 ).^(-0.5);
    aDOWN=( 1+dxDOWN.^2+dyDOWN.^2 ).^(-0.5);
    Unew=diffusion2(J,aLEFT,aDOWN,fmatrix,g);
    deltaU=max(max(abs(U-Unew)));    U=Unew;
    disp(['  iteration ' num2str(m) ',  max change in U = ' num2str(deltaU)])
    if deltaU<TOL, break, end
end

mesh(x,y,U')
xlabel x, ylabel y, hidden off, title('approximate minimal surface')