Lots of things to contemplate to get this to work. But the first step is realising your initial guess is currently fields of zeros yielding a singularity in the residual (and therefore the NaNs). So the following change will help with the singularity in the residual.
fns = Function(W)
fns.interpolate(Constant((1.0, 1.0, 1.0, 1.0)))
With regards to all the other joys of nonlinear problems maybe this will help.