I would like to solve the following system:

```
-∇ϕ = E
div(E) = ρ
0 = ∂ₜρ = div(μEρ) + s
```

It can be interpreted as the steady state of charged particle density ρ, generated by a process s

and moving around under under an external potential ϕ as well as its own electric field.

I tried the following:

```
from fenics import *
mesh = UnitIntervalMesh(5)
class ManufacturedSolution():
x = SpatialCoordinate(mesh)[0]
phi = x**2
E = -grad(phi)
rho = div(E)
mu = 2
s = div(E*mu*rho)
truth = ManufacturedSolution()
cell = mesh.ufl_cell()
deg = 1
fe_rho = FiniteElement("CG", cell, deg)
fe_E = VectorElement("CG", cell, deg)
fe_phi = FiniteElement("CG", cell, deg)
fe = MixedElement([fe_rho, fe_E, fe_phi])
V = FunctionSpace(mesh, fe)
v_rho, v_E, v_phi = TestFunctions(V)
u = Function(V)
rho, E, phi = split(u)
V_rho = V.sub(0)
V_E = V.sub(1)
V_phi = V.sub(2)
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0.0)
bc_rho = DirichletBC(V_rho, project(truth.rho, V_rho.collapse()), LeftBoundary())
bc_phi = DirichletBC(V_phi, project(truth.phi, V_phi.collapse()), "on_boundary")
bcs=[bc_rho, bc_phi]
F_rho = inner(truth.s - div(truth.mu*E*rho), v_rho) * dx
F_E = inner(div(E) - rho, v_phi)*dx
F_phi = inner(-grad(phi) - E, v_E) * dx
F = F_rho + F_phi + F_E
def relerror(phi_exact, phi):
scale = norm(phi_exact, "L2")
return errornorm(phi_exact, phi, "L2")/scale
solve(F == 0, u, bcs=bcs)
rho,E, phi = u.split()
print(relerror(project(truth.rho, V_rho.collapse()), rho))
print(relerror(project(truth.E , V_E .collapse()), E ))
print(relerror(project(truth.phi, V_phi.collapse()), phi))
plot(rho)
```

But I get

```
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 4.045e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
...
```

- I am looking for both being able to solve the issue for the specific example at hand and some general advice.
- Why does the solver not converge?
- Is the system instable? In the linear case, I would look at the assembled matrix or the condition number. Here I don’t really know how to start investigating the problem.
- How to reformulate the problem in a way that can be solved?