Based on the Navier-Stokes equations and a few parameterizations, the horizontal steady-state wind u(z) within a forest of height H satisfies:

$a(\frac{du}{dz})^2 + b\frac{du}{dz} \frac{d^2u}{dz^2} + cu +d\frac{du}{dz} + eu^2 + f= 0$, for $0<z<H$.

The coefficients a to f vary with the altitude z and are given initially (we can differentiable and integrate them as many times as needed).

At ground level: $u(z = 0) = 0$, $\frac{du}{dz}(z=0) = 0$.

At canopy top: $u(z = H) = U_H$, $\frac{du}{dz}(z=H) = K$ (constant)

I am trying to solve this equation for u(z) using a finite difference scheme, it would be great if someone could help me:

1/ Are Finite Differences even a good approach for this kind of problem ?

2/ If I rewrite the equation using the classical expressions $\frac{du}{dz} = \frac{u_{n+1}-u_{n-1}}{2h}, \frac{d^2u}{dz^2} =$ etc... I obtain square terms like $u_{i+1}u_{i-1}$ and I do not how what to do from there.

3/ I do not know how to use the Newton method or the Picard method correctly, is there a better way to rewrite the equation ? Using variables like $v = \frac{du}{dz}$ for example ?

At that point, I am not even sure if I am missing something obvious or if this is a really hard problem, any help would be greatly appreciated.

Thanks