I think one way to proceed is to use the Maclaurin expansions of $\sin$ and $\cos$: $\sin x \approx x - \frac{x^3}6$ and $\cos x \approx 1 - \frac{x^2}2$ when $x$ is small. If you plug in these approximate polynomials into your expressions, and then apply the sum-of-arguments formulas, polynomials will come out that will tell you the approximate values.
We would normally start by disregarding the second term of each polynomial, taking $\sin x \approx x$ and $\cos x \approx 1$, but if we do that here we get $\phi_1(x)\approx 0$ and $\phi_2(x)\approx 0$, which is correct, but less accurate than we want, and indicates that we disregarded too much.
The answer you linked to also links to an even older answer (with 100 upvotes!) which uses the same method. The earlier author also recommends a certain book citation; presumably if you really want more depth on this topic then that book would be a good source. Anyway, I'll refer to the answer you linked as "A1". I'll refer to the answer A1 links to as "A2".
I think A2's work is fully correct and well explained. A1's work is also basically correct, but they wrote the definition of $x_n$ backward (as you noted) and also made no effort to carefully track or bound the error terms (instead just removing them from the equation). Here I'll present a hopefully better explained version of A1's argument.
To start, we first define $x_n = u_n^{-2}$; note that's "$=$" not "$\approx$". Then the asymptotic recurrence becomes
$$x_{n+1} = 4 x_n + \frac 8 3 - \frac 4 {15} x_n^{-1} + O(x_n^{-2})$$
Following A1's lead, we drop the $- \frac{4}{15} x_n^{-1}$ term to make the recurrence easier to solve, getting:
$$
\begin{align}
x_{n+1} &= 4 x_n + \frac 8 3 + O(x_n^{-1}) \\
x_{n+1} + \frac 8 9 &= 4\left(x_n + \frac 8 9\right) + O(x_n^{-1}) \\
&= 4 \left( 4\left(x_{n-1} + \frac 8 9\right) + O(x_{n-1}^{-1}) \right) + O(x_{n}^{-1}) \\
&\, \, \,\vdots \\
x_{n} + \frac 8 9 &= 4^n \left(x_0 + \frac 8 9 \right) + O\left(\sum_{k=0}^{n-1} 4^{k} x_{n-1-k}^{-1}\right) \\
\end{align}
$$
To help simplify the error term, we can apply the above equation inductively to prove $x_n \sim 4^n x_0 $ as $x_0 \to \infty$. [EDIT: See bottom for a proof of this fact.] Note it makes sense to consider $x_0 \to \infty$ because we're aiming to end up with an asymptotic expansion valid as $u_0 \to 0$. The new estimate lets us write
$$
\begin{align}
O\left(\sum_{k=0}^{n-1} 4^{k} x_{n-1-k}^{-1}\right)
&= O\left(\sum_{k=0}^{n-1} 4^{k} (4^{n-1-k} x_0)^{-1}\right) \\
&= O \left( 4^{-n-1} x_0^{-1} \sum_{k=0}^{n-1} 4^{2k} \right) \\
&= O \left( 4^{-n} x_0^{-1} 16^n \right) \\
&= O \left( 4^{n} x_0^{-1} \right)
\end{align}
$$
Now plug this back into the main equation and replace $x_n \to u_n^{-2}$:
$$
\begin{align}
x_{n} + \frac 8 9 &= 4^n \left(x_0 + \frac 8 9 \right) + O \left( 4^{n} x_0^{-1} \right) \\
4^{-n} u_n^{-2} &= u_0^{-2} + \frac 8 9 - 4^{-n} \frac 8 9 + O \left( u_0^{2} \right)
\end{align}
$$
And finally, solving for $u_n$:
$$
\begin{align}
u_n &= \frac{2^{-n} u_0}{\sqrt{1 + \frac 8 9 (1-4^{-n}) u_0^2 + O \left( u_0^{4} \right)}}
\end{align}
$$
This finally matches the result A1 got, except that I've kept track of the error term so now we have an asymptotic estimate for how good this approximation is.
Takeaway
My overall takeaway is that A1 did know what they were doing, but they skipped a lot of steps esp involving error terms / set up the definition of $x_n$ backward / generally didn't bother explaining things very carefully. If I got more familiar with the technique then I could probably skip steps safely too, which would be nice because then I could confidently get the right result without needing so much writing. But for now I need all the steps because otherwise I might easily end up handwaving away some term that actually matters!
Final comment: A1's upgraded estimate
After this argument, A1 offers a second upgraded approximation. That calculation follows the same strategy as above, except without dropping the $- \frac 4 {15} u_n^2$ term from the original asymptotic recurrence relation. This makes the computations messier but ends up buying you one extra term in the final formula, which is why the error term ends up as $O(u_0^6)$.
If you wanted, you could keep computing more terms in the same way: Write the original asymptotic recurrence relation out a bit further, then apply this method to solve the asymptotic recurrence.
EDIT: Detailed proof of error bound lemma
Here are more details justifying this claim from above: "We can apply the above equation inductively to prove $x_n \sim 4^n x_0 $ as $x_0 \to \infty$."
The base case $n=0$ is trivial. Assuming the result for $k < n$, we have
$$
\begin{align}
x_{n} + \frac 8 9 &= 4^n \left(x_0 + \frac 8 9 \right) + O\left(\sum_{k=0}^{n-1} 4^{k} x_{n-1-k}^{-1}\right) \\
x_{n} &= 4^n x_0 + 4^n \frac 8 9 - \frac 8 9 + O\left(\sum_{k=0}^{n-1} 4^{k} 4^{k+1-n} x_0^{-1}\right) \\
\frac{x_{n}}{4^n x_0} &= 1 + \frac 8 9 x_0^{-1} - \frac 8 9 4^{-n} x_0^{-1} + 4^{-n} x_0^{-1} O\left(4^{1-n} x_0^{-1} \sum_{k=0}^{n-1} 4^{2k} \right) \\
\frac{x_{n}}{4^n x_0} &= 1 + \frac 8 {9x_0} - \frac 8 9 4^{-n} x_0^{-1} + x_0^{-2} O\left(1 \right) \\
\end{align}
$$
All RHS terms except for the $1$ will go to $0$ as $x_0 \to \infty$, uniformly in $n$ (we can pick one set of constants that work for all $n$), so $x^n \sim 4^n x_0$ as claimed.
Best Answer
The iteration has the form $$u_{n+1}=a_1u_n+a_3u_n^3+...$$ As usual in such situations (See the answer in Convergence of $\sqrt{n}x_{n}$ where $x_{n+1} = \sin(x_{n})$ with citation of de Bruijn: "Asymptotic Methods ..."), one can try a Bernoulli-like approach and examine the dynamics of $u_n^{-2}$. There one finds $$ \frac1{u_{n+1}^2}=\frac4{u_n^2(1-\frac13u_n^2+\frac15u_n^4\mp...)^2} =\frac4{u_n^2}+\frac83-\frac4{15}u_n^2+O(u_n^4)\tag1 $$ Thus for a first approximation use $$x_{n+1}=4x_n+\frac83\iff x_{n+1}+\frac89=4(x_n+\frac89)$$ so that $$u_n^{-2}\sim x_n=4^n(x_0+\frac89)-\frac89.\tag2$$
This gives as first approximation $$ u_n\sim \frac{2^{-n}u_0}{\sqrt{1+\frac89u_0^2(1-4^{-n})}}.\tag3 $$
For the next term use $v_n=(u_n^{-2}+\frac89)^{-1}$ and express (1) in terms of $v_n$ $$ \frac1{v_{n+1}}=\frac4{v_n}-\frac4{15}\frac{v_n}{1-\frac89v_n}+O(v_n^2) \tag4 $$ so that $$ \frac1{v_{n+1}}-\frac{4}{15^2}v_{n+1}=\frac4{v_n}-\frac1{15} v_n - \frac{1}{15^2}v_n+O(v_n^2)=4\left(\frac1{v_{n}}-\frac{4}{15^2}v_{n}\right)+O(v_n^2) \tag5 $$ and consequently $$ \frac1{v_{n}}-\frac{4}{15^2}v_{n}=4^n\left(\frac1{v_{0}}-\frac{4}{15^2}v_{0}+O(v_0^2)\right) \tag6 $$ As $\frac1v-\frac{4}{15^2}v=\frac1v(1-\frac4{15^2}v^2)$, leaving out the second term adds an error $O(v_n^2)$ which is a small fraction of $O(v_0^2)$. Thus $$ \frac1{u_n^2}+\frac89=\frac1{v_n}=4^n\left(\frac1{v_{0}}-\frac{4}{15^2}v_{0}+O(v_0^2)\right)=4^n\left(\frac1{u_0^2}+\frac89-\frac{4}{15^2}\frac{u_0^2}{1+\frac89u_0^2}+O(u_0^4)\right)\tag7 $$ so that the improved approximation is $$ u_n=\frac{2^{-n}u_0}{\sqrt{1+\frac89u_0^2(1-4^{-n})-\frac{4}{25}\frac{u_0^4}{9+8u_0^2}+O(u_0^6)}} \tag8 $$