Hi all,

I’d like to solve this system of PDE :

I ended up with this kind of formulation :

```
from fenics import *
mesh = UnitIntervalMesh(10)
V = VectorFunctionSpace(mesh, 'P', 1, 2)
u = Function(V)
u_n = Function(V)
v = TestFunction(V)
sol = list(split(u))
prev_sol = list(split(u_n))
test = list(split(v))
F = (sol[0] - prev_sol[0])/2*test[0]*dx - 1 * test[0]*dx
F += (sol[0]*sol[1] - prev_sol[0]*prev_sol[1])/2*test[1]*dx - 1 * test[1]*dx
solve(F == 0, u, [])
```

But this doesn’t converge (residual is nan). I reckon there’s something I’m not doing right here. Thanks for the help !

Rem