clear

Nrealiz=1000
kappa=0.1
dt=0.01
endt=10
Nt=floor(endt/dt)+1
histstart=-2;
histend=2;
histstep=0.1;
holdendx=zeros(Nrealiz,1);

tic
for i=1:Nrealiz
    x=zeros(Nt,1);
    t=zeros(Nt,1);
    x(1)=0; % initial condition
    t(1)=0;
    r=randn(Nt,1);
    for n=1:(Nt-1)
        x(n+1)=x(n)+(-x(n)^3+x(n))*dt+sqrt(2*kappa*dt)*r(n);
        t(n+1)=t(n)+dt;
    end
    holdendx(i)=x(Nt);
end
binedges=histstart:histstep:histend;
hist_unnorm=hist(holdendx,binedges);
pdf=hist_unnorm/(length(holdendx)*histstep);
figure
plot(binedges,pdf,'b.');

toc