////////////////////////////////////////////////////////////////////////////////
////        Equations de Dyson                                             ////
////       (champ moyen répulsif, avec schéma d'euler simple)              ////
///////////////////////////////////////////////////////////////////////////////



clf();
stacksize("max");
T=1;
n=100; //Pas de temps T/n
time_mesh=(T/n)*[0:n]';
N=100;  //Nombre de particules 
       //Exemples pour voir le demi-cercle : T=1, n=100 et N=1000 





//////////////////////////////////////////////////////////////
//Condition initiale (exemple) = la loi du demi-circle ///////
/////////////////////////////////////////////////////////////

U1=grand(1,N,'def');
U2=grand(1,N,'def');
X=gsort(2*sqrt(U1).*cos(2*%pi*U2),'lc','i');  


/////////////////////////////////////////////////////////////////


for t=1:n


Y=[];

for i=1:N
U=X(t,i)-X(t,:);
V=((U>0)|(U<0));
Y=[Y,(1/N)*sum(U(V)^(-1))];
end

////////////////////////////////////////////////////////////////////////////
////// COMPLETER la ligne rouge :  EVOLUTION DU SCHEMA d'EULER          ////
////////////////////////////////////////////////////////////////////////////

X=[X;	  ];

end



// AFFICHAGE DES RESULTATS


clf(0);xset("window",0);show_window();//Tracé des 100 premières trajectoires
    
    for i=1:100
        plot2d(time_mesh,X(:,i),i);
    end

clf(1);xset("window",2);show_window(); //Tracé de l'histogramme des états terminaux
histplot(linspace(-2*sqrt(T),2*sqrt(T),20),X(n,:),style=5);
x= linspace(-2*sqrt(T),2*sqrt(T))';
plot2d(x,(1/(2*%pi*T))*sqrt(4*T-x^2),style=3)