int n1 = 20;
mesh Th = square(n1,n1);
fespace Xh(Th,P1);
Xh u,v,delta,u0,w,uw;
int n = Xh.ndof;
real mu = 1e-2;
real rho = 1;
real[int] ret(n),aux(n),deltav(n);
u0 = x*(x-1)*y*(y-1); // initial guess
func real[int] F(real[int] & u)
{
  Xh uloc;
  uloc[] = u;
  varf Fvar(unused,v) = int2d(Th)(mu*(dx(uloc)*dx(v)+dy(uloc)*dy(v))
				  +rho*(uloc^2*v)
				  -v)
    +on(1,2,3,4,unused=0);
  ret = Fvar(0,Xh);
  return ret;
}
u = u0;
func real[int] JF(real[int] & delta)
{
  Xh deltaloc;
  deltaloc[] = delta;
  varf JFvar(unused,v) = int2d(Th)(mu*(dx(deltaloc)*dx(v)+dy(deltaloc)*dy(v))
				   +2*rho*(deltaloc*u*v))
    +on(1,2,3,4,unused=0);
  ret = JFvar(0,Xh);
  return ret;
}
real epsilon = 1e-4;
w = x^2*(x-1)^3*y^3*(y-1)^4;
cout << "J_F(u)w                       = " << JF(w[]).l2 << endl;
uw = u+epsilon*w;
aux = (F(uw[])-F(u[]));
aux /= epsilon;
cout << "(F(u+epsilon*w)-F(u))/epsilon = " << aux.l2 << endl;