%  
%this code determines the probability of decay vs. exit at x=L, starting
%from x=0
clear all
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

L = 1; % absorb at x=L
D = 1; % diffusion coefficient

Ng = 50;
for g = 1:Ng
gamma = 0.1*g; %decay rate
gm(g) = gamma;
dt = .001;  %time step
diff = sqrt(2*D*dt); % diffusion step size
 
N = 1000;  %number of particles
nd = 0;  % number of decayed particles
for j = 1:N % for each particle
    
    x = 0;
   decay = 0;
        
    while ((abs(x)<L)& (decay ==0)) % step until x reaches L or particle decays
        x = abs( x+diff*randn(1,1)); % move by diffusion  (abs give reflection at x=0)
        
        if(gamma*dt>rand(1,1))
            decay = 1; %this particle has decayed
            nd = nd +1; % count the number that decay
        end
        
    end
end

 pd(g) =  1-nd/N ;  % the percentage that reached x=L
end

plot(gm,pd,'*',gm, 1./cosh(sqrt(gm)),'--','linewidth',2)
xlabel('\alpha L^2/D','fontsize',20)
ylabel('Probability of reaching x = L','fontsize',20)