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\Big(\frac{du}{dz}\Big)^{\!2} + b\frac{du}{dz} \frac{d^2u}{dz^2} + cu +d\frac{du}{dz} + eu^2 + f= 0\:\text{ 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, \quad\frac{du}{dz}\Big|_{z=0} = 0.
$$
At canopy top:
$$
u|_{z = H} = U_H, \quad\frac{du}{dz}\Big|_{z=H} = K\text{ (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:

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

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.

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