/////////////////////////////
// Agglomération de globules rouges
//////////////////////////////
//
//Paramètres-discrétisation
K=15000; h=0.0002; T=K*h; gama=0.5; landa=10; a=0.004, N=50
t=[0:h:T]
//
//Fonction phi
deff('expre=phi(d)','expre=[(gama/d)*log(d/a)]')
//
//Génération aléatoire du vecteur initial
x=rand(1,1,1)
X0=[]
for i=1:N
    x=2*rand(1,1,1)-1
    X0=[X0;(i+0.1*x)/(N+1)]
end
X1=X0
//
//Matrice d'approximation
Mapp=[X0,X1]

//La fonction PHI
function y = PHI(x)
y=[]
y=[y;phi(x(2)-x(1))-phi(x(1))] 
for i=2:N-1
    y=[y;phi(x(i+1)-x(i))-phi(x(i)-x(i-1))] 
end
y=[y;phi(1-x(N))-phi(x(N)-x(N-1))] 
end

//Le schéma numérique
for i=1:K
    Xi=2*Mapp(:,i+1)-Mapp(:,i)+h*h*PHI(Mapp(:,i+1))-landa*h*(Mapp(:,i+1)-Mapp(:,i))
    Mapp=[Mapp,Xi]
end

t=[t,T+h]

//Tracé
figure(1);
clf()
for i=1:N
    plot(t,Mapp(i,:),'b-')
end
xlabel('Temps (s)'); 
ylabel('Position');


