function [U,uex]=heatbm(x0,t,dt,N)
% HEATBM  Solution of the time dependent heat equation
%    u_t = (1/2) u_xx
% on entire real line by simulation of Brownian motion.  Initial condition
% is         /  1,   0 < x <1,
%     f(x) = |
%            \  0,   otherwise.
% Format:
%    [U,uex] = heatbm(x0,t,dt,N)
% where
%      U   = approx solution to heat equation
%      uex = exact solution to heat equation
%      x0  = points x at which to compute u(x,t); also locations of start-
%               ing points of Brownian particles;  x0 may be a row vector
%      dt  = time step
%      N   = number of Brownian particles for simulation; use N>=2
% ELB 5/10/05
% See also:  BROWN,  BROWNPOT

x=repmat(x0,N,1);  J=ceil(t/dt);  dt=t/J;
dx=sqrt(dt);  % scaling for Brownian motion
for j=1:J,  x = x + dx*randn(size(x));  end
U=sum(f(x))/N;  % compute average;  U has same size as x0
uex=exact(x0,t);  % compute exact solution

function z=f(x);   z=(x>0)&(x<1);

function u=exact(x,t);
u=0.5*erfc((x-1)/sqrt(2*t))-0.5*erfc(x/sqrt(2*t));