%solving the Kolmogorov forward eq.

dt=0.0001;
av=1; %dx^2/dt=av
dx=sqrt(av*dt);

T=15; %final time
L=5; %space interval is [-L, L]
a=-L; b=L; 
Nt=floor(T/dt);
Nx=floor((b-a)/dx); 
ceni=floor(-a/dx);
xv=[0:Nx-1];
xv=a+dx*xv;

rho=zeros(1,Nx);
rhon=rho;
rhof=rho; rhofr=rho; rhofl=rho;
rhod=rho; rhodr=rho; rhodl=rho;

isig=0.05;
imu=0;
rho=1/(isig*sqrt(2*pi))*exp(-(xv-imu).^2/(2*isig^2)); %approx. to delta f.
maxinit=max(rho);

plot(xv, rho, 'Color',[1 0 0]);

stop=0; 
i=1;
while (stop==0)
 
    rhof=(ff1(xv).^2).*rho; 
    rhofr(1:Nx-1)=rhof(2:Nx); rhofr(Nx)=rhof(1);
    rhofl(2:Nx)=rhof(1:Nx-1); rhofl(1)=rhof(Nx);
    rhod=fd1(xv).*rho; 
    rhodr(1:Nx-1)=rhod(2:Nx); rhodr(Nx)=rhod(1);
    rhodl(2:Nx)=rhod(1:Nx-1); rhodl(1)=rhod(Nx);
    rhon=rho-dx*(rhodr-rhodl)/(2*av)+0.5*(rhofr+rhofl-2*rhof)/av;
     
    rho=rhon;
   
    if (mod(i, 100)==0) %draw
        plot(xv, rho, 'Color',[1 0 0], 'LineWidth', 2);
        str=sprintf('t=%2.4f', i*dt);
        text(-L+1,1,str,...
          'HorizontalAlignment','left');
        set(gca, 'YLim',[0 2]);
        set(gca, 'XLim',[-L L]);
        set(gca, 'LineWidth',2);
        drawnow;
    end
    i=i+1;
    if ((i>Nt)||(max(rho)>2*maxinit)) 
        stop=1;
    end
end
