function [U,err]=brownpot(n,x0,y0,rho,Deltax);
% BROWNPOT  Solves potential equation problem
%    u_xx + u_yy = 0   on  D=(0,1)x(0,1),  u=g(x,y) on bdry D
% at (x0,y0) by *simulated Brownian motion* on n particles 
% starting at (x0,y0).  Simulation is by random walk with 
% normally distributed steps of size Deltax; default for Deltax 
% is 0.01.   Compares to exact solution
%    u(x,y) = x^3 - 3 x y^2; 
% the same function gives boundary values g(x,y).
% Form:
%    [U,err]=brownpot(n,x0,y0,rho,Deltax)
% where  (last two are input arguments are optional):
%         U = approximation to u(x0,y0)
%       err = U - u(x0,y0)
%         n = number of particles in simulation
%   (x0,y0) = coordinates of starting point
%      [rho = fraction which need to hit; 0 <= rho <= 1; default 0.5]
%   [Deltax = distance of step of random walk; default 0.01]
% Example 1:
%    >> [U,err]=brownpot(1000,0.5,0.5)
% Example 2:
%    >> [U,err]=brownpot(1000,0.5,0.5,0.99)
% Example 3:
%    >> [U,err]=brownpot(1000,0.5,0.5,0.99,0.1)
% ELB 4/11/03

if nargin==3  % use defaults depending on number of input arguments
  rhon=0.5*n;  dx=0.01;
elseif nargin==4, 
  rhon=rho*n;  dx=0.01;
elseif nargin==5, 
  rhon=rho*n;  dx=Deltax;
else, error('needs 3, 4 or 5 input arguments'), end
maxsteps = ceil(1/dx^2);

x = repmat(x0,n,1);  y = repmat(y0,n,1);  % all start at (x0,y0)
hit = false(n,1);  % none have hit boundary yet
numhits = 0;  sumg = 0; tic
for k=1:maxsteps
   newhit = ~hit & ( (x<=0) | (x>=1) | (y<=0) | (y>=1) );
   if any(newhit)
      x(newhit&(x<=0))=0;  x(newhit&(x>=1))=1;  % move back to bdry
      y(newhit&(y<=0))=0;  y(newhit&(y>=1))=1;
      sumg = sumg + sum(g(x(newhit),y(newhit)));  % evaluate g at new point
      numhits = numhits + length(x(newhit));
      if numhits >= rhon, break, end; %stop when rho fraction have hit bdry
      hit = hit | newhit;  end
   x = x + dx*randn(n,1);  y = y + dx*randn(n,1);  % new positions
   if toc>1, fprintf(1,'.'), tic, end
end
fprintf(1,'\n')
if k>=maxsteps, warning('exceeded max number of steps'), end
U=sumg/numhits;
err=U-g(x0,y0);

function z=g(x,y);  z = x.^3 - 3 * x .* y.^2;  %  bdry values