border a(t=0,2*pi) {x = cos(t); y = sin(t); label=1;}
mesh disk = buildmesh(a(50));

plot(disk, wait=true);

fespace Vh(disk, P1);
Vh u, v;

solve pb(u,v) = 
      int2d(disk)(dx(u)*dx(v)+dy(u)*dy(v))
     +int2d(disk)(-v)
     +on(1, u=0);
plot(u, value=1, wait=1);
Vh err=u-(1-x^2-y^2)/4;
plot(err, value=1, wait=1);