Solve functional equation $ h(y)+h^{-1}(y)=2y+y^2 $ I was trying to solve a certain physics problem, and encountered the functional equation that contains a function $h$ and its inverse $h^{-1}$:
\begin{equation}
h(y)+h^{-1}(y)=2y+y^2.\tag{1}
\end{equation}

Q: Does $(1)$ have a unique solution and is it possible to find it in closed form?

Equation $(1)$ looks quite simple and probably it is simple to analyze too, but I couldn't figured out how.

Background information: The physics problem I was trying to solve is finding the dependence of the current $J_0$ on the voltage $U_0$ in this infinite chain that contains nonlinear elements with quadratic volt-ampere characteristics $I(V)=\alpha V^2$ and ohmic resistors $R$:

According to dimensional analysis one can write
$$
J_0(U_0)=\frac{1}{\alpha R^2}f(\alpha RU_0).
$$
Solving simple system of equations I obtained the following functional equation for the unknown function $f$
$$
f(x)=(x-f(x))^2+f(x-f(x)).\tag{2}
$$
Now introduce another function $h$ according to
$$
x-f(x)=h(x).
$$
Then equation $(2)$ becomes
$$
x-h(x)=h^2(x)+h(x)-h(h(x)).
$$
Let $h(x)=y$, then $x=h^{-1}(y)$ and we get $(1)$.
This is not a textbook problem and I don't even know whether it has a solution. I was studying it out of curiosity.
 A: From
$$
x = 2h(x)+h^2(x)-h(h(x))\Rightarrow h(x) = \pm\sqrt{h(h(x))+x+1}-1
$$
now using an iterative approximation over the positive leaf such as
$$
h_{k+1}(x) = \sqrt{h_k(h_k(x))+x+1}-1
$$
or according to the MATHEMATICA script which follows, (inspired after a fruitful discussion with Semiclassical)
Clear[h]
h[x_, 1] := x
h[x_, n_] := h[x, n] = Sqrt[h[h[x, n - 1], n - 1] + x + 1] - 1

we obtain beginning with $h_1(x) = x$ the successive approximations in red, and black for $h_5(x)$ showing a good convergence for $x > 0$. 

Also the accomplishment for $h^{-1}(x)+h(x) = 2x+x^2$ can be analyzed in the following plot
err5 = y - h[2 y + y^2 - h[y, 5], 5];
Plot[Abs[err5], {y, 0, 4}, PlotStyle -> {Thick, Black}]


A: The following is not intended as an answer so much as it is a development of the physics problem. (Hence why it's a wiki answer.) Let $U_k$ be the potential difference across the $k$th nonlinear element with corresponding current $\alpha U_k^2$, for $k\geq 1$. By virtue of Kirchoff's loop rule, we deduce that the voltage across the $k$th resistor is given by $U_{k-1}-U_{k}$. Since each resistor is ohmic, we conclude that the current in the $k$th resistor is $$J_{k-1}=R^{-1}(U_{k-1}-U_{k}).$$ In particular, we have  $J_0=R^{-1}(U_{0}-U_{1})$. Kirchoff's junction rule then demands $$J_{k} = J_{k-1}-\alpha U_k^2.$$
Introducing $(u_k,j_k):=(\alpha R U_k,\alpha R^2 J_k)$ and rearranging, we obtain the dimensionless recurrence relations 
\begin{align}
u_k = u_{k-1}-j_{k-1},\qquad
j_k &= j_{k-1}-u_k^2\\&=j_{k-1}-(u_{k-1}-j_{k-1})^2.
\end{align}
This may be compactly expressed as $$(u_k,j_k)=g(u_{k-1},j_{k-1})=\cdots=g^k(u_0,j_0)$$ where $g(u,j):=(u-j,j-(u-j)^2)$.  That is, the sequence of voltages and currents is obtained by iterating from an initial choice of $(u_0,j_0)$.
However, on physical grounds we are only concerned with solutions for which the currents and voltages are positive and monotonically decreasing to zero. This is rather delicate, as numerical experimentation demonstrates that the fixed point $(u,j)=(0,0)$ is badly unstable. (I don't know how to formally prove this. We can note, though, that the conditions $u-j>0$ and $j>j-(u-j)^2>0$ require $j<u<j+\sqrt{j}$. Hence the range of possible $j$ gets smaller and smaller as $j\to 0^+$, which to me seems consistent with the origin being unstable.)
Finally, to obtain the advertised functional equation, we seek a solution of the form $j_k=f(u_k)$. In that case, applying the first equation to $j_{k+1}=f(u_{k+1})$ yields
$$j_{k+1} = f(u_{k+1}) = f(u_k-j_k) = f(u_k-f(u_k)).$$
On the other hand, from the second equation we have
$$j_{k+1} = j_k - u_{k+1}^2 = j_k - (u_k-j_k)^2 = f(u_k)-(u_k - f(u_k)^2.$$
Together, we obtain the desired functional equation $$f(x-f(x)) = f(x)-f(x-f(x))^2$$ under the identification $x=u_k$. (From this we may further conclude that $h(u_k) = u_{k+1}$, i.e., the voltages are obtained by iterating $h$ on $u_0$.)
I'm not confident I know how to justify the condition $j_k=f(u_k)$ rigorously. Physically, however, it has a simple interpretation: Suppose we draw a vertical line between the first nonlinear element and the second resistor. Then on the right we have a copy of the infinite chain, but subject to voltage $V_1$ and $J_1$. Since it's the same chain, it must have the same voltage-current relationship as the original, i.e., $j_1=f(u_1)$. We may then repeat this logic with the next element-resistor pair and so on, yielding $j_j=f(u_k)$ for all $k$ as desired.
At this point I altogether run out of firm conclusions. But I do have some more observations:


*

*In a problem with multiple dimensionful parameters, it's often wise to study limiting cases for which the problem simplifies. For instance, one has the trivial limit $U_k,J_k\to 0$ as $U_0\to 0$. More interesting is the case $\alpha\to \infty$, where each nonlinear element is a short circuit and so the entire current $J_0$ will flow through the first nonlinear element (path of least resistance). Therefore $J_0=R^{-1} U_0$ as $\alpha R U_0\to\infty$. Other limits are not helpful: If $\alpha\to 0$, then all the nonlinear elements become open circuits and therefore one has an infinite chain of identical resistors, i.e., infinite resistance. Similarly, if $R\to 0$ then the voltage across each nonlinear element is $U_0$ and therefore the current required would be infinite. Neither situation is physical and so there's no evident conclusion to draw. (To render them physical, one could either introduce ohmic resistances on the branches or consider a finite chain.)

*It's worth noting that, while the dimensionless variables chosen above seem obvious enough, they're not the only ones possible. For instance, we could just have well have taken $u'_k:=U_k/U_0$, $j'_k:=(R/U_0)J_k$ to obtain the dimensionless equations
$$j'_{k-1} = u'_{k-1} - u'_k,\qquad j'_k = j'_{k-1} - \gamma {u'_k}^2$$
where $\gamma = \alpha R U_0$. The principal benefit of this is that we can now explicitly consider $\alpha \to \infty$, in which case $\gamma\to \infty$ and so the second equation collapses to $u'_k=0$ for $k\geq 1$. Hence $j'_1=u_0$ as stated previously. This moreover raises the possibility of solving the equations perturbatively in powers of $\gamma^{-1}$, though I've run into trouble proceeding along this line.

*In the above, I've presented the problem in terms of a system of nonlinear first-order difference equations. This is readily converted into a single nonlinear difference equation of second-order:
$$u_k^2 = j_{k-1}-j_k = (u_{k-1}-u_k)-(u_{k}-u_{k+1})=u_{k+1}-2u_k+u_{k-1}.$$
This second-order difference equation is analogous to the differential equation $u(x)^2 = 2u''(x)$ with $j(x):=-u'(x)$. By taking a first integral, we obtain $$\frac{1}{3}u(x)^3 = u'(x)^2+C.$$  The requirement that $u(x),j(x)\to 0$ as $x\to \infty$ then imposes $C=0$. The resulting first-order ODE is separable with solution $u(x)=u(0)(1+x\sqrt{u(0)/12})^{-2}$. Then $j(0)=u(0)^{3/2}/\sqrt{3}$, suggesting $$j_0\propto \frac{u_0^{3/2}}{\sqrt{3}}\implies J_0\propto \frac{(\alpha R U_0)^{3/2}}{\sqrt{3}\alpha R^2}$$ (for small $U_0$, I think?) as the asymptotics for the original difference equation.  The above is not exactly rigorous, so I don't know how much stock to put in it; however, it does seem to match up with what numerics I've done.

*One solution idea, in line with the above use of a first integral to solve the differential equation, is to search for a conservation law, i.e., a function $H(x,y)$ such that $$H(u_k,j_k) = (H\circ g)(u_{k-1},j_{k-1})=H(u_{k-1},j_{k-1}).$$ (Or, for $u_k$ alone, a function $H'$ such that $H'(u_{k+1},u_k)=H'(u_k,u_{k-1})$.) This would dictate $$H(u_0,j_0)=\cdots =H(u_k,j_k)=\cdots=\lim_{k\to \infty} H(u_k,j_k)=H(0,0).$$ As such,  the set of physical $(u_0,j_0)$ would be prescribed as the level set of $H$ through $(0,0)$. However, I've yet to come up with such $H(x,y)$ and so have no definitive conclusions here.
A: I tried to find a more general chain that has an exact solution and then obtain the chain in the question as a limiting case of this more general problem. However, it turned out that for exact solution of this more general problem there needs to be a constraint for parameters. Because of this constraint this approach didn't work out.
Consider the infinite network in the picture:

Here $I_X(V)=\alpha\sqrt{V}$ and  $I_Y(V)=\beta V^2$ are two nonlinear elements.  One can combine the constants $R$, $\alpha$, $\beta$ to obtain a quantity that has dimension of voltage
\begin{equation}
    V_0=(\alpha/\beta)^{2/3},
\end{equation}
and a dimensionless parameter
\begin{equation}
    \varepsilon=R(\alpha^2\beta)^{1/3}.
\end{equation}
Then
\begin{equation}
    \alpha=\varepsilon\sqrt{V_0}/R,\quad \beta=\varepsilon/(V_0R).
\end{equation}
To simplify calculations let's introduce shorthand $Z$ for the part of the network as shown in picture a) below.

One can easily calculate its equivalent volt-ampere characteristic
\begin{equation}\label{Z}
    I_Z(V)=\frac{V_0}{R}F_Z\left(\frac{V}{V_0}\right),\quad \text{где} \quad F_Z(t)=\varepsilon t^2+\varepsilon\sqrt{\varepsilon^2+t}-\varepsilon^2.
\end{equation}
According to dimensional analysis we can write for the VAC of the infinite chain
\begin{equation}\label{AB}
    I_{AB}(U_0)=\frac{V_0}{R}F\left(\frac{U_0}{V_0}\right),
\end{equation}
where $F$ - is an unknown dimensioneless function of a dimensionless argument. We want to find this function.
We find the following equation for $F$:
\begin{equation}\label{f}
    F(t)=F(t-F(t))+F_Z(t-F(t)).
\end{equation}We want to find a solution $0<F(t)<t$.
Introducing 
\begin{equation}\label{h}
    H(t)=t-F(t),
\end{equation}
we get
\begin{equation}\label{hfinal}
    H(t)+H^{-1}(t)=2t+F_Z(t).\tag{*}
\end{equation}
Assume the solution is of the form
$$
F(t)=t+\frac{R}{V_0}(a-\sqrt{a^2+bV_0t})
$$
with some constants $a$ and $b$.
Then $$
H(t)=\frac{R}{V_0}(\sqrt{a^2+bV_0t}-a),
$$
and
$$
H^{-1}(t)=\frac{V_0t^2+2aRt}{bR^2}.
$$ 
Substituting these into $(*)$ we get
$$
\frac{R}{V_0}(\sqrt{a^2+bV_0t}-a)+\frac{V_0t^2+2aRt}{bR^2}=2t+\varepsilon t^2+\varepsilon\sqrt{\varepsilon^2+t}-\varepsilon^2.\tag{**}
$$
One can easily see that in order for this equation to satisfy identically for all $t>0$ one needs to satisfy the conditions
$$
\varepsilon=1,\quad a=V_0/R, \quad b=V_0/R^2.
$$
Unfortunately, replacing the resistance $2R$ with more general $r$ doesn't lead to anything new, because the requirement analogous to $(**)$ can only be satisfied when $r=2R$.
Thus we get that we can find the VAC of this infinite chain under the constraint $\alpha^2\beta R^3=1$:
$$
I_{AB}(U_0)=\frac{U_0+V_0-\sqrt{V_0^2+V_0U_0}}{R}, \quad \text{where} \quad V_0=(\alpha/\beta)^{2/3}.
$$
Unfortunately, the requirement 
$$
\alpha^2\beta R^3=1
$$
spoils everything and it is not possible to use this solution to find the solution of the initial question as a limiting case.
