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

cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

 %  splitting probability
 N = 1000; % number of particles
 K = 19;  % number of starting positions
 L = 1; % size of the domain
D = 1; % diffusion coefficient
dt = .0005;  % time step
diff = sqrt(2*D*dt);
 dx = L/(K+1);
 for jx = 1:K
     xs = jx*dx;  % starting positions
     xst(jx) = xs;
 for j = 1:N
    clear x X
    x = xs;
    k = 1;
    
    X(1) = 0;
    while ((x<L)&(x>0))  % step until x = 0 or L
        x =  (x+diff*randn(1,1));
        k = k + 1;
        X(k) = x; 
        
    end
    Xs(j) = x; % record the final location
    
end
    q(jx) = sum((Xs>=L))/N;  % calculate the percentage that exited at L
 end
 x = [0:0.1:L];
 piL= x/L; %This is the theoretical answer
 
    figure(1)
    plot(xst,q,'*',x,piL,'--','linewidth',2)
    
    xlabel('starting position','fontsize',20)
    ylabel('\pi_L','fontsize',24)
    
     if(cmykflg==1)
print('../../figs_c/chapt_4/split_prob','-deps','-cmyk')
  end

 
    al = sum((q-0.5).*(xst-L/2));
    bt = sum((xst-L/2).^2);
    a = al/bt;
    figure(2)
    plot(xst,q,'*',x,piL,'k--',x,a*(x-L/2)+0.5,'r--','linewidth',2)
    
    xlabel('starting position','fontsize',20)
    ylabel('\pi_L','fontsize',20)