N=100;
h=1/N;
M=10000;
Z=randn(M,N);
X=zeros(M,N+1);

a=0.2;
sigma=0:0.1:5;
p=zeros(size(a));

for j=1:length(sigma)
   for i=1:N
      X(:,i+1)=X(:,i)+a*h*X(:,i).*(2-X(:,i))+sigma(j)*sqrt(h)*Z(:,i);
   end
   p(j)=sum(max(X,[],2)<1)/M;
end
plot(sigma,p)
xlabel('sigma')
ylabel('p')