mesh Th=square(100,100);
 fespace Vh(Th,P1);     // P1 FE space
 Vh uh,vh,u1=0;              // unkown and test function. 
 func f=1;                 //  right hand side function 
 func g=0;                 //  boundary condition function
 real dt =0.01;
 int i=0;
 problem Poisson(uh,vh,init=i) =                    //  definion of  the problem 
    int2d(Th)( uh*vh+dt*(dx(uh)*dx(vh) + dy(uh)*dy(vh)) ) //  bilinear form
  - int2d(Th)( (u1+dt*f)*vh  )                          //  linear form
  + on(1,2,3,4,uh=g) ;                      //  boundary condition form

real[int] colorhsv=[  // color hsv model
  4./6., 1 , 0.5, // dark blue
  4./6., 1 , 1, //  blue
  5./6., 1 , 1, //  magenta
  1, 1. , 1, //  red
  1, 0.5 , 1 // light red 
   ];

  real[int]  viso(20);
  for (int i=0;i<viso.n;i++)
  viso[i]=i*0.005;

for (i=0;i<10;i++)
 {
  Poisson; // solve the problem 
  plot(uh,value=true,hsv=colorhsv,viso=viso,fill=1);  // to see the result
  u1=uh; 
}