% In this code we follow diffusing/decaying particles to determine the
% distribution of decay positions
clear

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

L=1;
mu = 1 % decay rate, dimensionless

d = 2 ; % a parameter
dt = 0.01;  % small time step
N = 0;
Ntot=1000;
for k = 1:Ntot
    d = .05*k;
    sc = sqrt(2*d*dt);

    dk(k) = d;
    N = 0;
for j = 1:1000
x = 0; % initialize the position
 
while(abs(x)<1)
r = rand(1);
if r>dt;
    x = x +sc*randn(1);
else
    N = N + 1;
   x = 2;
end
 
end

end
Nk(k) = N/Ntot;
N
end 

f = 1-1./cosh(1./sqrt(dk));

plot(dk,Nk,dk,f,'--') 