int n=40, m=n/5;
real eps=0.001;

real x0=0, x1=5;
real y0=0, y1=1; 
mesh Th=square(n,m, [x0+(x1-x0)*x, y0+(y1-y0)*y]);
plot(Th, wait=1);


fespace Vh(Th, P1);
Vh u,v;
Vh betax=4*y*(1-y);

solve convdiff(u,v) =
     int2d(Th)(eps*dx(u)*dx(v)+eps*dy(u)*dy(v)+betax*dx(u)*v)
     +on(4, u=0)
     +on(1,3,u=1);
plot(u);