Unfortunately, your proposed functional equation:
$$ h(y)+h^{-1}(y)=2y+y^2 \tag{1} $$ has no
differentiable solutions near $\,y=0.\,$ Suppose that
$$ h(y) = 0 + a_1 y + O(y^2). \tag{2} $$
Then its inverse function is
$$ h^{-1}(y) = 0 + \frac1{a_1} y + O(y^2). \tag{3} $$
Substituting equations $(2)$ and $(3)$ into $(1)$ gives
$$ h(y)\!+\!h^{-1}(y) = \frac{a_1^2+1}{a_1}y \!+\! O(y^2)
= 2y \!+\! O(y^2). \tag{4}$$
The only solutions for $\,a_1\,$ are $\,a_1=1\,$ and
$\,a_1=-1.\,$ However, we are given that
$$ 0<h(x)<x \tag{5}$$
and thus the only possible value is $\,a_1=1.\,$
Suppose we want more terms in the power series
$$ h(y) = 0 + y + a_2y^2 + O(y^3). \tag{6}$$
The inverse function is now
$$ h^{-1}(y) = 0 + y - a_2y^2 + O(y^3). \tag{7}$$
Adding these two equations together is the equation
$$ h(y) + h^{-1}(y) = 0 + 2y + 0y^2 + O(y^3). \tag{8}$$
This implies that equation $(1)$ can not be satisfied.
There is a possibility that there exists an exponent
$\,e\,$ not an integer such that
$\, h(y) = 0 + y + c y^e + \cdots \,$ which perhaps
should be studied according to comments.
In fact, define $$ g(x) : = \sqrt{h(x^2)}. \tag{9}$$
Then solving for a power series expansion gives
$$ g(x) = x - \frac{x^2}{\sqrt{6}} + \frac{x^3}6 + O(x^4) \tag{10}$$ which implies that
$$ h(x^2) = (x^2+x^4/2)+\frac1{\sqrt{6}}f(x) \tag{11} $$
where
$$ f(x) \!=\! -2 x^3
\!-\!\frac{11x^5}{24} \!+\! \frac{117x^7}{1280}
\!-\!\frac{5491x^9}{110592} \!+\! \frac{156538363x^{11}}{3715891200} \!+\! O(x)^{13}. \tag{12}$$
A Wolfram Language code to calculate $\,g(x)\,$ is
ClearAll[x, g, gx];
gx[3] = x - 1/Sqrt[6]*x^2 + O[x]^3;
Do[g = Normal[gx[n]] + O[x]^(2+n); gx[n+1] = Simplify[
g + (x^2 + (Normal[g]/.x -> g)^2 - g^4 - 2*g^2) * 3 /
((4+n)*x^2*Sqrt[6])],
{n, 3, 6}]
As some comments suggest, it seems that the power series
for $\,f(x)\,$ has zero radius of convergence. This makes
finding properties of it difficult. Perhaps we need to let
$\,x\,$ approach $\,\infty.$ In this case we find that
$\,h(x) \approx \sqrt{x}\,$ with infinite number of other
terms. Using equation $(1)$ we can find the expansion
$$ h(x) = x^{1/2} - 1 + \frac12 x^{-1/4} + \frac1x s(x) \tag{13} $$
where $$ s(x) :=
\sum_{k=2}^\infty 2^{-k}(x^{-\frac32 2^{-k}} -
x^{-2^{-k}} ). \tag{14} $$
This gives moderately good approximations as $\,x\,$ gets large and down to $1.$
One method that leads to equations $(13)$ and $(14)$
is as follows. From equation $(1)$ we immediately get
$$ h^{-1}(y) = y^2 + 2y - h(y) \tag{15} $$
and if $\,y=h(x),\,$ then
$$ x = h(x)^2 + 2h(x) - h(h(x)). \tag{16} $$
We start with an approximation and try to find what
additional term will satisfy equation $(16)$. So we guess
$$ h(x) = x^{1/2} + cx^e + \cdots \tag{17} $$
where $\,\cdots\,$ denotes terms with smaller exponents.
Substitute equation $(17)$ in equation $(16)$ to get
$$ x \!=\! (x \!+\! 2cx^{1/2+e} \!+\! \cdots) \!+\!
2x^{1/2} \!+\! \cdots \!=\! x \!+\! 2x^{1/2}(cx^e \!+\! 1)
+ \cdots. \tag{18} $$
This implies $\,0 = cx^e + 1.\,$ Solving this we get
our next guess as
$$ \,h(x) = x^{1/2} - 1 + cx^e + \cdots. \tag{19} $$
Repeating this process leads to equations $(13)$ and
$(14)$.
The series in $(14)$ appears
to converge but I have no proof, only numerical evidence.
I also have no proof that $(13)$ satisfies the functional
equation $(1)$.
An answer to this question by 'Semiclasical' contains
sequence recursions
$$ u_n = u_{n-1} - j_{n-1} \quad \text{ and } \quad
j_n = j_{n-1} - u_n^2. \tag{20} $$
with the property that $\, u_n = h(u_{n-1})\,$ and
that $\,u_n \to 0^+.\,$ I found that for suitable
starting values of $\,u_1\,$ and $\,j_1\,$ we have
$$ u_n = 6n^{-2} - \frac{15}2n^{-4} + \frac{663}{40}
n^{-6} - \frac{43647}{800}n^{-8} + \cdots. \tag{21} $$
An answer to this question by 'Cesareo' uses recursion to construct a sequence of functions $\,h_n(x)\,$ which seems to converge to a global solution. It would be nice to give a proof of this.
Best Answer
I would like to give a solution of both the instances of the problem which appeared here in different times, i.e. the version posted before the OP edited his post and the current one. First I want to solve the problem as it is stated by the OP.
Problem I
$$f(y+f(x))=f(x)f(y)+f(f(x))+f(y)-xy.$$ As already noticed, setting $y=0$, one arrives to $$f(f(x))=f(0)f(x)+f(f(x))+f(0)$$ which leads to $f(x)\equiv -1$, which is not a solution of the problem, or $f(0)=0$. Incidentally, even though it is not crucial for the solution, notice that if we assume that there is an $\bar x$ such that $f(\bar x)=0$, setting $x=\bar x$ we infer that $f(y)=f(y)-\bar xy$. In particular, setting $y=1$, we conclude $\bar x=0$.
Next, plug inthe equation $f(y)$ instead of $y$ to get $$f(f(y)+f(x))=f(x)f(f(y))+f(f(x))+f(f(y))-xf(y).$$ Interchange the roles of $x$ and $y$ to obtain similarly $$f(f(x)+f(y))=f(f(x))f(y)+f(f(y))+f(f(x))-yf(x).$$ Compare these two equations to conclude that $$\frac{f(f(x))+x}{f(x)}=\frac{f(f(y))+y}{f(y)}=k,$$ where $k$ is a suitable constant.
Substitute in the original equation to get $$f(y+f(x))=f(x)f(y)-x+kf(x)+f(y)-xy.$$ Denote $a:=f(-1)$ and choose $y=-1$ to obtain $$f(f(x)-1)=(a+k)f(x)+a$$
Edit Thanks to Dejan who pointed out a flaw in my proof
Plug in the original equation $f(y)-1$ instead of $y$ to obtain $$\begin{aligned}f(f(y)+f(x)-1)&=f(x)[(a+k)f(y)+a]-x+kf(x)-x(f(y)-1)\\&=(a+k)(f(x)+f(y))-xf(y).\end{aligned}$$ Exchange $x$ and $y$ to obtain $$f(f(y)+f(x)-1)=(a+k)(f(x)+f(y))-yf(x).$$ Compare these two equations to deduce that $$\frac{f(x)}{x}=\frac{f(y)}{y}=a.$$
Substitute in the original equation to find out that the admissible values for $a$ are just $\pm 1$ i.e. $$f(x)=\pm x.$$
Problem II
Now let's see what happens if the equation is written instead as $$f(x+f(y))=f(x)f(y)+f(f(x))+f(y)-xy.$$ Setting $x=y=0$ one deduces that $$f(0)^2+f(0)=0$$ from which $f(0)\in \{-1, 0\}$.
Assume then $f(0)=0.$ Then plug in the equation $x=0$ to deduce that $f(y)=f(f(y)).$
Let's see that $f$ is injective. Suppose there are $y\neq z$ with $f(y)=f(z)$. Use the equation with $y$ and $z$ and deduce that in this case $x(y-z)=0$ should hold for all $x\in\mathbb R$, which is plainly impossible. Then $f$ is injective. Now, recall that in our hypothesis the equation states as $$f(x+f(y))=f(x)f(y)+f(x)+f(y)-xy.$$ Exchange $x$ and $y$ to obtain $$f(y+f(x))=f(y)f(x)+f(y)+f(x)-yx.$$ It follows from the injectivity of $f$ that $$f(x)-x=f(y)-y=k=0.$$ It is easy to verify that the identity is indeed a solution of our equation.
Assume now $f(0)=-1$. Then, choosing $x=0$, one deduces $$f(f(y))=f(f(-1))=a.$$ Substitute then $x$ with $f(x)$ to derive $$f(f(x)+f(y))=af(y)+f(a)+f(y)-f(x)y.$$ Similarly, also $$f(f(x)+f(y))=af(x)+f(a)+f(x)-f(y)x$$ holds. Therefore we finally find that $$\frac{f(x)}{(a+1+x)}=\frac{f(y)}{(a+1+y)}=c.$$ It follows that there are constants $c,d\in\mathbb R$ such that $$f(x)=cx+d.$$ From the fact that $f(f(x))=a$, we deduce $c=0$, and from $f(0)=-1$, it follows $f\equiv -1$ which is not a solution of the problem.