int n1 = 20;
mesh Th = square(n1,n1);
fespace Xh(Th,P1);
Xh u,v,delta,u0,w;
int n = Xh.ndof;
real mu = 1e-2;
real rho = 1;
real tol = 1e-12;
real[int] ret(n),rhs(n),deltav(n);
u0 = x*(x-1)*y*(y-1); // initial guess
u = u0;
// first method
problem AR(delta,v) = int2d(Th)(mu*(dx(delta)*dx(v)+dy(delta)*dy(v))
			      +2*rho*(delta*u*v))
  +int2d(Th)(mu*(dx(u)*dx(v)+dy(u)*dy(v))
	     +rho*(u^2*v)
	     -v)
  +on(1,2,3,4,delta=0);
AR;
cout << "delta " << delta[].linfty << endl;
while (delta[].linfty > tol)
  {
    u = u+delta;
    AR;
    cout << "delta " << delta[].linfty << endl;
  }
cout << "norm(u) = " << u[].l2 << endl;
// second method
u = u0;
varf Fvar(unused,v) = int2d(Th)(-mu*(dx(u)*dx(v)+dy(u)*dy(v))
				-rho*(u^2*v)
				+v)
  +on(1,2,3,4,unused=0); // notice the sign
ret = Fvar(0,Xh);
varf JFvar(delta,v) = int2d(Th)(mu*(dx(delta)*dx(v)+dy(delta)*dy(v))
				+2*rho*(delta*u*v))
  +on(1,2,3,4,delta=0);
matrix A = JFvar(Xh,Xh,solver=CG); // here you can specify the solver
delta[] = A^(-1)*ret;
while (delta[].linfty > tol)
  {
    u = u+delta;
    varf Fvar(unused,v) = int2d(Th)(-mu*(dx(u)*dx(v)+dy(u)*dy(v))
				    -rho*(u^2*v)
				    +v)
      +on(1,2,3,4,unused=0); // notice the sign
    ret = Fvar(0,Xh);
    varf JFvar(delta,v) = int2d(Th)(mu*(dx(delta)*dx(v)+dy(delta)*dy(v))
				    +2*rho*(delta*u*v))
      +on(1,2,3,4,delta=0);
    matrix A = JFvar(Xh,Xh,solver=CG); // here you can specify the solver
    delta[] = A^(-1)*ret;
    cout << "delta " << delta[].linfty << endl;
  }
cout << "norm(u) = " << u[].l2 << endl;
// third method
u = u0;
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;
}
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;
}
rhs = F(u[]);
deltav = 0; // initial guess
LinearCG(JF,deltav,rhs,eps=-1.e-30); // ||rhs-JF*deltav||_2^2<-eps
cout << "deltav " << deltav.linfty << endl;
while (deltav.linfty > tol)
  {
    u[] = u[]+deltav;
    rhs = F(u[]);
    LinearCG(JF,deltav,rhs,eps=-1.e-30); // ||rhs-JF*deltav||_2^2<-eps
    cout << "deltav " << deltav.linfty << endl;
  }
cout << "norm(u) = " << u[].l2 << endl;