 mesh Th=square(10,10);
 fespace Vh(Th,P1);     // P1 FE space
 Vh uh,vh,u;              // unkown and test function. 
 func f=1;                 //  right hand side function 
 func g=0.;                 //  boundary condition function
 
 problem laplace(uh,vh,solver=GMRES,tgv=1e5) =                    //  definion of  the problem 
  - int2d(Th)( -dx(uh)*dx(vh) - dy(uh)*dy(vh) ) //  bilinear form
  - int1d(Th,1)( -uh*vh) 
  - int1d(Th,1)(  vh)
  - int2d(Th)( f*vh )                          //  linear form
  + on(2,3,4,uh=g) ;                      //  boundary condition form

 problem laplacep(uh,vh,solver=CG,tgv=1e5) =                    //  definion of  the problem 
    int2d(Th)( dx(uh)*dx(vh) + dy(uh)*dy(vh) ) //  bilinear form
  + int1d(Th,1)( uh*vh) 
  + int1d(Th,1)( - vh)
  + int2d(Th)( -f*vh )                          //  linear form
  + on(2,3,4,uh=g) ;                      //  boundary condition form

  
  laplace; // solve the problem plot(uh); // to see the result
  u=uh;
  laplacep; // solve the problem plot(uh); // to see the result
  
  plot(uh,u,value=true);
  u[] -= uh[];
  cout << "Diff min = "<< u[].min << " max=" << u[].max << endl ;
