Why the plot function doesn't plot the velocity?


the plot([uite1, uite2], pite, coef=0.2, cmm=" [Ux, Uy] - p", WindowIndex=1);
only plot the result of the first step, then it only plot the contour of the region, how to fix this?, the full code is

load "iovtk"
verbosity=0;
// Parameters
int nn = 15;
int nnPlus = 5;
real l = 1.;
real L = 15.;
real hSlope = 0.1;
real H = 6.;
real h = 0.5;

real beta = 0.01;

real eps = 9.81/303;
real nu = 1.;
real numu = nu/sqrt(0.09);
real nuep = pow(nu,1.5)/4.1;
real dt = 0.;
// Mesh
border b1(t=0, l){x=t; y=0;}
border b2(t=0, L-l){x=1.+t; y=-hSlope*t;}
border b3(t=-hSlope*(L-l), H){x=L; y=t;}
border b4(t=L, 0){x=t; y=H;}
border b5(t=H, h){x=0; y=t;}
border b6(t=h, 0){x=0; y=t;}

mesh Th=buildmesh(b1(nnPlus*nn*l) + b2(nn*sqrt((L-l)^2+(hSlope*(L-l))^2)) + b3(nn*(H + hSlope*(L-l))) + b4(nn*L) + b5(nn*(H-h)) + b6(nnPlus*nn*h));
plot(Th);

// Fespaces
fespace Vh2(Th, P1b);
Vh2 Ux, Uy,uite1,uite2,Uxst,Uyst;
Vh2 Vx, Vy,v1,v2;
Vh2 Upx, Upy=0;
fespace Vh(Th,P1);
Vh p=0, q,pp=p,pite,vp,pst;
Vh Tp, T=35,Tite,vt,Tst;
Vh k=log(0.27), kp=k,kite,vk,kst;
Vh ep=log(0.002104), epp=ep,epite,vep,epst;
 
fespace V0h(Th,P0);
 V0h muT=1;
 V0h prodk, prode;
 Vh kappa=0.25e-4, stress;
// Macro
 macro grad(u) [dx(u), dy(u)] //
 macro Grad(u1,u2) [dx(u1),dy(u1), dx(u2),dy(u2)] //
 macro Gradt(u1,u2) [dx(u1),dx(u2), dy(u1),dy(u2)] //
 macro Div(u1,u2) (dx(u1) + dy(u2)) //
macro UgradV(u1,u2,v1,v2) [[u1,u2]'*[dx(v1),dy(v1)],[u1,u2]'*[dx(v2),dy(v2)]]//
macro ugradv(u1,u2,v) [u1,u2]'*grad(v)//
// Problem
real alpha = 0.;


// Initialization with stationary solution


// Initialization
dt = 0.05; 
Tp = T - 10*((x<1)*(y<0.5) + (x>=1)*(y+0.1*(x-1)<0.5));
Upx=0+(x==0)*(y>=0)*(y<=h)*3;
plot(Upx,  cmm=" Ux ", wait=1);
plot(Upy,  cmm=" Ux ", wait=1);
// Convergence loop
real T0 = clock();
cout  << " - dt = " << dt << endl;
   alpha = 1/dt;

   // Time loop
   real t = 0.;
   for (int i = 0; i <= 500; i++){
      uite1=Upx;
      uite2=Upy;
      pite=pp;
      Tite=Tp;
      kite=kp;
      epite=epp;
      t=t+dt;
      cout << "Time step " << i << " - t = " << t << endl;
      for (int n=0;n<=4;n++){
         solve xstar ([Uxst,Uyst,pst,Tst,kst,epst],[v1,v2,vp,vt,vk,vep])
         =int2d(Th,optimize=0)(alpha*([uite1,uite2]'*[v1,v2]-Tite*vt-kite*vk-epite*vep)
         +(UgradV(uite1,uite2,uite1,uite2))'*[v1,v2]
         +(UgradV(uite1,uite2,uite1,uite2))'*[v1,v2]
         +ugradv(v1,v2,pite)
         +(0.09)*exp(2*kite-epite)*(2*kite-epite)*(Grad(uite1,uite2)+(Gradt(uite1,uite2)))'*Grad(v1,v2)
         +(0.09)*exp(2*kite-epite)*(Grad(uite1,uite2)+Gradt(uite1,uite2))'*Grad(v1,v2)
         +(9.81/303)*Tite*v2-Div(uite1,uite2)*vp-vk*(ugradv(uite1,uite2,kite)+ugradv(uite1,uite2,kite))-exp(epite-kite)*(epite-kite)*vk+(0.09)*exp(2*kite-epite)*(2*kite-epite)*[dx(kite),dy(kite)]'*[dx(kite),dy(kite)]*vk
         +(0.09)*exp(2*kite-epite)*2*[dx(kite),dy(kite)]'*[dx(kite),dy(kite)]*vk
         -(0.09)*exp(2*kite-epite)*(2*kite-epite)*[dx(vk),dy(vk)]'*[dx(kite),dy(kite)]-(0.09)*exp(2*kite-epite)*[dx(vk),dy(vk)]'*[dx(kite),dy(kite)]
         +(0.09)*exp(kite-epite)*(4*dx(uite1)*dx(uite1)+4*dy(uite2)*dy(uite2)+2*(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vk
         +(0.09)*exp(kite-epite)*(kite-epite)*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vk
         -vep*(ugradv(uite1,uite2,epite)+ugradv(uite1,uite2,epite))
         -(1.92)*exp(epite-kite)*(epite-kite)*vep
         +(0.09/1.3)*exp(2*kite-epite)*(2*kite-epite)*[dx(epite),dy(epite)]'*[dx(epite),dy(epite)]*vep
         +(0.09/1.3)*exp(2*kite-epite)*2*[dx(epite),dy(epite)]'*[dx(epite),dy(epite)]*vep
         -(0.09/1.3)*exp(2*kite-epite)*(2*kite-epite)*[dx(epite),dy(epite)]'*[dx(vep),dy(vep)]
         -(0.09/1.3)*exp(2*kite-epite)*[dx(epite),dy(epite)]'*[dx(vep),dy(vep)]
         +((1.44*0.09))*exp(kite-epite)*(4*dx(uite1)*dx(uite1)+4*dy(uite2)*dy(uite2)+2*(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vep
         +((1.44*0.09))*exp(kite-epite)*(kite-epite)*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vep
         -vt*(ugradv(uite1,uite2,kite)+ugradv(uite1,uite2,kite))
         -0.1*exp(2*kite-epite)*(2*kite-epite)*[dx(vt),dy(vt)]'*[dx(Tite),dy(Tite)]
         -0.1*exp(2*kite-epite)*[dx(vt),dy(vt)]'*[dx(Tite),dy(Tite)]
         )
         +int1d(Th,b1,b2,b5,optimize=0)(
            (0.09/1.3)*2.4952*exp(0.5*kite)*(0.5*kite)*vep
         )+int1d(Th,b1,b2,optimize=0)(
            0.41*((0.09)^0.25)*(exp(0.5*kite)*(0.5*kite)*uite1*v1+exp(0.5*kite)*uite1*v1)
         )-int2d(Th,optimize=0)(
            alpha*([uite1-Upx,uite2-Upy]'*[v1,v2]-(Tite-Tp)*vt-(kite-kp)*vk-(epite-epp)*vep)
            +(UgradV(uite1,uite2,uite1,uite2))'*[v1,v2]
            +(grad(pite))'*[v1,v2]
            +0.09*exp(2*kite-epite)*(Grad(uite1,uite2)+(Gradt(uite1,uite2)))'*Grad(v1,v2)
            +(9.81/303)*(Tite-35)*v2-(Div(uite1,uite2))*vp
            -(ugradv(uite1,uite2,kite))*vk
            -exp(epite-kite)*vk
            +(0.09)*vk*(exp(2*kite-epite)*grad(kite)'*grad(kite))
            -(0.09)*(exp(2*kite-epite)*grad(kite)'*grad(vk))
            +(0.09)*exp(kite-epite)*vk*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))
            -(ugradv(uite1,uite2,epite))*vep
            -1.92*exp(epite-kite)*vep
            +(0.09/1.3)*vep*(exp(2*kite-epite)*grad(epite)'*grad(epite))
            -(0.09/1.3)*(exp(2*kite-epite)*grad(epite)'*grad(vep))
            +(1.44*0.09)*exp(kite-epite)*vep*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))
            -(ugradv(uite1,uite2,Tite))*vt
            -0.1*(exp(2*kite-epite)*grad(Tite)'*grad(vt))
         )-int1d(Th,b1,b2,b5, optimize=0)(
            (0.09/2.6)*2.4952*exp(0.5*kite)*vep+(0.09/2.6)*2.4952*exp(0.5*kp)*vep)
         -int1d(Th,b1,b2, optimize=0)(0.41*((0.09)^0.25)*exp(0.5*kite)*uite1*v1)
         -int2d(Th, optimize=0)(alpha*([Uxst,Uyst]'*[v1,v2]-Tst*vt-kst*vk-epst*vep)
         +(UgradV(Uxst,Uyst,uite1,uite2))'*[v1,v2]
         +(UgradV(uite1,uite2,Uxst,Uyst))'*[v1,v2]
         +ugradv(v1,v2,pst)
         +(0.09)*exp(2*kite-epite)*(2*kst-epst)*(Grad(uite1,uite2)+(Gradt(uite1,uite2)))'*Grad(v1,v2)
         +(0.09)*exp(2*kite-epite)*(Grad(Uxst,Uyst)+Gradt(Uxst,Uyst))'*Grad(v1,v2)
         +(9.81/303)*Tst*v2-Div(Uxst,Uyst)*vp-vk*(ugradv(Uxst,Uyst,kite)+ugradv(uite1,uite2,kst))-exp(epite-kite)*(epst-kst)*vk+(0.09)*exp(2*kite-epite)*(2*kst-epst)*[dx(kite),dy(kite)]'*[dx(kite),dy(kite)]*vk
         +(0.09)*exp(2*kite-epite)*2*[dx(kite),dy(kite)]'*[dx(kst),dy(kst)]*vk
         -(0.09)*exp(2*kite-epite)*(2*kst-epst)*[dx(vk),dy(vk)]'*[dx(kite),dy(kite)]-(0.09)*exp(2*kite-epite)*[dx(vk),dy(vk)]'*[dx(kst),dy(kst)]
         +(0.09)*exp(kite-epite)*(4*dx(uite1)*dx(Uxst)+4*dy(uite2)*dy(Uyst)+2*(dy(uite1)+dx(uite2))*(dy(Uxst)+dx(Uyst)))*vk
         +(0.09)*exp(kite-epite)*(kst-epst)*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vk
         -vep*(ugradv(Uxst,Uyst,epite)+ugradv(uite1,uite2,epst))
         -(1.92)*exp(epite-kite)*(epst-kst)*vep
         +(0.09/1.3)*exp(2*kite-epite)*(2*kst-epst)*[dx(epite),dy(epite)]'*[dx(epite),dy(epite)]*vep
         +(0.09/1.3)*exp(2*kite-epite)*2*[dx(epite),dy(epite)]'*[dx(epst),dy(epst)]*vep
         -(0.09/1.3)*exp(2*kite-epite)*(2*kst-epst)*[dx(epite),dy(epite)]'*[dx(vep),dy(vep)]
         -(0.09/1.3)*exp(2*kite-epite)*[dx(epst),dy(epst)]'*[dx(vep),dy(vep)]
         +((1.44*0.09))*exp(kite-epite)*(4*dx(uite1)*dx(Uxst)+4*dy(uite2)*dy(Uyst)+2*(dy(uite1)+dx(uite2))*(dy(Uxst)+dx(Uyst)))*vep
         +((1.44*0.09))*exp(kite-epite)*(kst-epst)*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vep
         -vt*(ugradv(Uxst,Uyst,kite)+ugradv(uite1,uite2,kst))
         -0.1*exp(2*kite-epite)*(2*kst-epst)*[dx(vt),dy(vt)]'*[dx(Tite),dy(Tite)]
         -0.1*exp(2*kite-epite)*[dx(vt),dy(vt)]'*[dx(Tst),dy(Tst)]
         )
         -int1d(Th,b1,b2,b5, optimize=0)(
            (0.09/1.3)*2.4952*exp(0.5*kite)*(0.5*kst)*vep
         )-int1d(Th,b1,b2, optimize=0)(
            0.41*((0.09)^0.25)*(exp(0.5*kite)*(0.5*kst)*uite1*v1+exp(0.5*kite)*Uxst*v1)
         )
         +on(b5,Uxst=0,Uyst=0)
         +on(b6,Uxst=3,Uyst=0,kst=log(0.27),epst=log(0.002104),Tst=25)+on(b2,Uyst=-Upx*N.x/N.y,Tst=30)+on(b1,Uyst=0,Tst=30);
         uite1=Uxst;
         uite2=Uyst;
         pite=pst;
         Tite=Tst;
         kite=kst;
         epite=epst;
      }
      Upx=uite1;
      Upy=uite2;
      pp=pite;
      Tp=Tite;
      kp=kite;
      epp=epite;
      plot([uite1, uite2], pite, coef=0.2, cmm=" [Ux, Uy] - p", WindowIndex=1);


   }
cout << "Total Time = " << clock()-T0 << endl;

Does the problem occur with a simpler (to follow) model.
I guess it is best if we include a simple example to illustrate the issue, rather than an example with lots of code lines which may not add much info to reproduce the issue.
Regards

This works, but it is an extremely inefficient implementation. If you coarsen up the mesh a bit, you will see the plots you want.