Error while trying to solve coupled PDEs

Keep in mind that the python layer of DOLFIN works with UFL so that many objects have an algebraic representation by inheriting the functionality of UFL. So when you solve the nonlinear problem, the Gateaux derivative is computed using UFL with the UFL representation of the DOLFIN objects (i.e. your Function).

Regarding convergence: The original error you were getting has nothing to do with convergence. I just noticed that running the code with fix above and the residual and Jacobian as defined for your formulation will require more care/work, since you have fractional powers.