Finite difference for a highly nonlinear equation - The wind within the forest 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
 A: As noted already in the comments, your boundary conditions seem off. Note that generically for a second-order BVP one expects to impose only two boundary conditions; you have 4. 
Once you’re sure you’ve formulated the problem correctly, finite differences are a good first approach.  A finite difference discretization will leave you with a system of nonlinear algebraic equations to solve, and I would suggest using a library for that (e.g. fsolve in Matlab or scipy).  You could also directly use a library that solves BVPs.
To get a basic understanding of this kind of thing, I recommend reading the first two chapters of LeVeque's book on finite differences which includes an example of what can happen if you impose boundary conditions that make the problem ill-posed. 
Since your problem is probably convection-dominated, you may find it necessary to use upwinded (one-sided) finite differences for the convective terms in order to avoid oscillations.
Finally, I recommend posting questions like this on scicomp.SE, as there are many experts that read that site but aren't on Mathoverflow.
