Good closed form approximation for iterates of $x^2+(1-x^2)x$
We have that $f(x) - x = x^2 - x^3$, so if we pretend that $n$ is a continuous variable, and treat the iterates $x(n)$ as functions of this continuous variable, then we are led to the differential equation:
$$\frac{dx}{dn} = x^2 - x^3$$
Which has the solution:
$$\log\left[\frac{x(n)}{1-x(n)}\right]-\frac{1}{x(n)} = \log\left[\frac{x(0)}{1-x(0)}\right]-\frac{1}{x(0)} + n$$
You can then approximately solve for $x(n)$, what is then clear is that for $x(n)$ near zero the $\frac{1}{x(n)}$ term dominates while for $x(n)$ near $1 $ the logarithmic term dominates.
FYI, here is a picture illustrating Count's answer at https://math.stackexchange.com/a/1872168/241. Blue curves are the zeroth through tenth iterates of $f$; red curves are obtained by numerically inverting $F(f^{\circ N}(x)) \approx N + F(x)$ for $F(x) := \log \frac{x}{1-x} + \frac{1}{x}$. I will mention that this is an improvement on the exact fourth-order approximation $f^{\circ N}(x) \approx x+Nx^2+(N^2-2N)x^3+([N-1]^3-2[N-1]^2-3[N-1])x^4$ near $0$, and I imagine an improvement on all exact fixed-order approximations.
On a separate note, I will also mention the Carleman approach I originally mentioned can be pursued along the lines of http://www.sciencedirect.com/science/article/pii/S0022247X98959868.
In blue: $g(x,25)$. In orange: $0.5\tanh(x(1+x)^{25})+0.5$. Your function does stay at $y=1$ thereafter.
The behaviour near zero is entirely influenced by $u(x)=x(1+x)$, neglecting the $-x^3$ term. The obvious approximate composition that I mentioned above, $u^n(x)\approx x(1+x)^n$ simply doesn't increase fast enough. Perhaps your Carlemann matrix aproach could be used on this simpler function to yield a more meaningful approximation near zero?
The $\tanh$ motivation should also be obvious. I hope this helps.