%% LANGEVIN integrates the Langevin equation
% using the standard Euler-Maruyama method. 
% The solution is Ornstein-Uhlenbeck noise.  
% Solution can be seen using plot(t,x);
% Histogram of solution can be seen using plot(rhox,rhoy);
%
%   Copyright 2003
%   Centre for Nonlinear Dynamics in Physiology and Medicine

tot=20000;  %total number of integration time steps per realizations
navgs=10;   %total number of realizations
dt=0.01;    % integration time step
trans=1000; %number of time steps to discard as transients
alf=1.0;    % alpha
dnz=1.0;    %intensity of the Gaussian white noise process; dnz=0.5*sigma^2
xinit=1.0;  %initial condition x(0)
t=ones(1,tot);  %time vector
bins=200;       %number of bins to use for histogram of solution 
rhoy=zeros(1,bins);  %histogram of the solution
transtime=trans+1;
sig=sqrt(2*dnz);
fac=sig*sqrt(dt);
for realiznumber=1:navgs
    realiznumber
    rng(sum(100*clock)); % initialize the random number generator
    t=zeros(1,tot);
    x=xinit*ones(1,tot); %initialize solution vector
    for i=2:tot
        t(1,i)=(i-1)*dt;
        x(1,i)=x(1,i-1)-alf*x(1,i-1)*dt+fac*randn;
    end
    rhoy=rhoy+hist(x(transtime:tot),bins);
end
[rhodum,rhox]=hist(x(transtime:tot),bins);