Hi everyone, thanks for reading!

I am trying to solve the following system of ODEs:

\left[ \left( \nu+ \sigma^* \frac{k}{\omega} \right) k_z \right]_z = \ \beta^* k \omega - \frac{k}{\omega} \frac{(z P + 1)^2}{(\nu + k/ \omega)^2} \\ \left[ \left( \nu + \sigma \frac{k}{\omega} \right) \omega_z \right]_z = \ \beta \omega^2 - \alpha \frac{(z P + 1)^2}{(\nu + k/ \omega)^2}

where everything apart from k and \omega are known constants, and k_z means dk/dz. I have Dirichlet boundary conditions k(0) = k(1) = 0 and \omega(0) = \omega(1) = \omega_{wall}

If I expand this and define dk = k_z and d \omega = \omega_z, I get four equations

dk = k_z \\
k_{zz} = \left( \beta^*k \omega - \sigma^*\frac{dk \cdot \omega - k \cdot d\omega}{\omega^2} dk - \frac{k}{\omega} \left( \frac{zP + 1}{\nu + k/\omega} \right)^2 \right) / (\nu + \sigma^* k/ \omega) \\
d \omega = \omega_z \\
\omega_{zz} = \left( \beta \omega^2 - \sigma \frac{dk \cdot \omega - k \cdot d\omega}{\omega^2} d\omega - \alpha \left( \frac{zP + 1}{\nu + k/\omega} \right)^2 \right) / (\nu + \sigma k/ \omega) \\

I am confident this is correct because it works when I solve it using scipy’s solve_bvp.

Below is my implementation in FEniCS. The Newton solver goes to nan in the 0th iteration:

`Newton iteration 0: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)`

Could someone please help me debug this? Thank you!

```
from dolfin import *
nu = Constant(1/58.4)
alpha = Constant(13/25)
beta = Constant(0.0708)
beta_star = Constant(9/100)
sigma = Constant(1/2)
sigma_star = Constant(3/5)
P = Constant(-2.001556)
ncells = 10000
d = 1/ncells
om_wall = Constant(6*nu/(beta*d**2))
mesh = IntervalMesh(ncells, 0, 1)
P2 = FiniteElement('Lagrange', interval, 2)
P1 = FiniteElement('Lagrange', interval, 1)
Welm = MixedElement([P2, P1, P2, P1])
W = FunctionSpace(mesh, Welm)
bcsys = [DirichletBC(W.sub(0), Constant(0.0), 'near(x[0], 0)'), # k = 0 at z = 0
DirichletBC(W.sub(0), Constant(0.0), 'near(x[0], 1)'), # k = 0 at z = 1
DirichletBC(W.sub(2), om_wall, 'near(x[0], 0)'), # \omega = \omega_{wall} at z = 0
DirichletBC(W.sub(2), om_wall, 'near(x[0], 1)')] # \omega = \omega_{wall} at z = 0
fns = Function(W)
k, kz, om, omz = split(fns) # k, k_z, \omega, \omega_z
q, r, s, t = split(TestFunction(W))
expr = Expression('x[0]', degree=1) # z
# dk = k_z
k_eqn = k.dx(0) - kz
# k_zz = ...
kz_eqn = kz.dx(0) - (
beta_star * k * om
- k/om / (nu + k/om)**2 * (expr*P + 1)**2
- sigma_star * (kz * om - k * omz) / om**2 * kz
) / (nu + sigma_star * k/om)
# d \omega = \omega_z
om_eqn = om.dx(0) - omz
# \omega_zz = ...
omz_eqn = omz.dx(0) - (
beta * om**2
- alpha / (nu + k/om)**2 * (expr*P + 1)**2
- sigma * (kz * om - k * omz) / om**2 * omz
) / (nu + sigma * k/om)
F = (k_eqn * q + kz_eqn * r + om_eqn * s + omz_eqn * t) * dx
Jac = derivative(F, fns, TrialFunction(W))
solve(F == 0, fns, J=Jac, bcs=bcsys)
```

Edit: added kz.dx(0) and omz.dx(0) terms in kz_eqn and omz_eqn, respectively.