Besoin de comprendre la logique de cet algo

macro model 2
macro typeFE() P1//
//macro DBC true// by default NBC 

verbosity=0; 
bool delay=true; //true si tau n'est pas nul
//load "UMFPACK64"

// Parameters
// k2=a, k3=beta, sigma2=sigma, sigma1=b
real L=2.;
real T=50.;
real a=10.;
real b=0.;
real D=1.e-3;
real beta=5.;
real sigma=0.4;
int Nbx=1e4;
real dt=.01;
real tau=10.;
real u0=1.;
real i0=0.;
real v0=10.;
real x0=0.1;
int Nbt=1e2;


real Dx=L/Nbx;

load "msh3"
meshL Th=segment(Nbx,[x*L,0.]);
fespace Vh2(Th,typeFE);
Vh2 Id = (x<=x0) ? v0 : 0.;

IFMACRO(model,2)
   Vh2 V,I,U,VH,IH,UH,V0=Id,I0=i0,U0=u0;
   macro SYSTEM
   {
      solve systemV(V,VH,solver=LU)=
         int1d(Th)(V*VH/dt + D*dx(V)*dx(VH) + sigma*V*VH)
         - int1d(Th)(V0*VH/dt + beta*Itmtau*VH)
      IFMACRO(DBC)
         + on(1,V=v0)
      ENDIFMACRO
         ;
      solve systemU(U,UH,solver=LU)=
         int1d(Th)(U*UH/dt + a*U*V*UH)
         - int1d(Th)(U0*UH/dt);
      solve systemI(I,IH,solver=LU)=
         int1d(Th)(I*IH/dt + b*I*IH)
         - int1d(Th)(I0*IH/dt + a*U*V*IH);
      V0=V; I0=I; U0=U;
   }//
ENDIFMACRO
Vh2 Itmtau=I0,Itausave;
int ntau=tau/dt;
real[int,int] Itau(Vh2.ndof,ntau+2);
Itausave=I0;
Itau(:,0)=Itausave[];
int op=0,ittau=0,ittmtau=0;

int it,itsave;
int Thnv = Th.nv;
fespace Vh(Th,P1);
Vh Vs=x,Is=x,Us=x;
real[int,int] Usave(Thnv,3),Vsave(Thnv,3),Isave(Thnv,3);
Usave(:,itsave)=Us[];
Isave(:,itsave)=Is[];
Vsave(:,itsave)=Vs[];
for (real t=0;t<T;t+=dt){
   if(t<tau){
      if(delay){
         Itmtau=0;
      } else {
         Itmtau=I0;
      }
      SYSTEM;
      if(delay){
         ittau++;//l'indice de la ligne d'après dans Itau
         Itausave=I;
         Itau(:,ittau)=Itausave[];
      }
   } else {
      if(delay){
         if(ittau==(ntau+1)){
            ittau=0;
            ittmtau=0;
         }
         Itmtau[] = Itau(:,ittmtau);
         SYSTEM;
         ittmtau++;
         Itausave=I;
         Itau(:,ittau)=Itausave[];
         ittau++;
      } else {
         Itmtau=I0;
         SYSTEM;
      }
   }
  
   if(!(it%max(1.,1./dt))){
      Vs=V;
      Is=I;
      Us=U;
      itsave++;
      Usave(:,itsave)=Us[];
      Isave(:,itsave)=Is[];
      Vsave(:,itsave)=Vs[];
      Usave.resize(Thnv,itsave+3);
      Isave.resize(Thnv,itsave+3);
      Vsave.resize(Thnv,itsave+3);
   }
   if(!(it%(max(1.,1./dt/2.)))){
      real Vmax,Vmin,Umax,Umin,Imax,Imin;
      Vs=V;
      Vmax=Vs[].max; Vmin=Vs[].min;
      Vs[]=Vs[]/Vmax;
      Is=I;
      Imax=Is[].max; Imin=Is[].min;
      Is[]=Is[]/Imax;
      Us=U;
      Umax=Us[].max; Umin=Us[].min;
      Us[]=Us[]/Umax;
      meshL ThmvV=movemesh(Th,[x/L,y+Vs]);//pour normaliser on met x/L (le domaine est entre 0 et 1), x/l*2 entre 0 et 2.
      meshL ThmvI=movemesh(Th,[x/L,y+Is]);
      meshL ThmvU=movemesh(Th,[x/L,y+Us]);
      plot(ThmvV,ThmvI,ThmvU,dim=2,wait=0,cmm="time: "+t+", U min = "+Umin+", U max = "+Umax+", V min = "+Vmin+", V max = "+Vmax+", I min = "+Imin+", I max = "+Imax);
   }
   op=1;
   it++;
}

ofstream foutU("u.txt");
ofstream foutI("i.txt");
ofstream foutV("v.txt");
for(int m=0;m<Th.nv;m++) {
   for(int n=0;n<itsave;n++) {
      foutV << Vsave(m,n) << " ";
      foutI << Isave(m,n) << " ";
      foutU << Usave(m,n) << " ";
   }
   foutV<<endl;
   foutI<<endl;
   foutU<<endl;
}