systems of equations – Numerically find the half-iterate of a quadratic


Let’s say we have a function $f(x)=n+mx+lx^2$, and want to find another function $h(x) = a+bx+cx^2 $ such that $h(h(x)) = f(x)$. Here’s what you get when composing h on itself:

$$h(h(x)) = a+ab+b^{2}x+bcx^{2}+a^{2}c+2abcx+2ac^{2}x^{2}+b^{2}cx^{2}+2bc^{2}x^{3}+c^{3}x^{4}$$

If we ignore the $x^3$ and onwards terms and gather the constants, $x$ terms, and $x^2$ terms, we can get a system of three equations which ignores enough terms to not be overly specified, but including enough terms to be a half-decent approximation:

$$a+ab+a^{2}c = n,$$
$$b^{2} + 2abc = m,$$
$$bc+2ac^{2}+b^{2}c = l$$

My original plan was to get $c$ isolated in each of these like so:
$$c = frac{(n – a (b + 1))}{a^2},$$
$$c = frac{(m – b^2)}{2 a b},$$
$$c = frac{l}{(b^2 + b)}$$

This way I could treat each of these as 3D functions and find their intersection (much like how you can solve a system of 2 equations by treating the equations as 2D functions). Using the identity that $2x-y-z=0$ if $x=y=z$, I could reduce the intersection problem further down into one equation (which can also be treated as a 3D function) which I need to find the zero/root of:
$$2cdot{frac{(n – a (b + 1))}{a^2}} -frac{(m – b^2)}{2 a b} – frac{l}{(b^2 + b)} = 0$$

The problem is, I don’t know of any root-finding algorithms supporting 2-argument functions, since for example, Netwon’s method requires a derivative which we can’t do here. So, is my approach wrong, is there a root-finding method I don’t know of, or is there no solution?