Based on the Navier-Stokes equations and a few parameterizations, the horizontal steady-state wind at a given altitude 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 (they are "well behaved", we can differentiate or integrate them as many times as needed).

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/ Is Finite Diff 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