No answer, just a comment to illustrate the functionality of the Schröder-function which is not familiar for several participiants here plus a (futile) attempt to a real-to-real solution (but possibly reflecting parts of the Kneser-method)
#1 - On the Schröder-mechanism
I applied the Schröder-mechanism which runs into a power series involving complex terms. Here is a graph for iterates in steps of 1/60 from a couple of starting points in the interval $0 \ldots \pi/2$ Near that two borders, the Schröder-function is difficult to handle and the iterates in that neighbourhoods are questionable.
For starting points below $\pi/4$ :
For starting points above $\pi/4$ :
Of course, the spirals/trajectories can be continued infinitely towards the fixpoint.
#2 - On an attempt for a real-to-real solution
I investigated the possibility to find a real-to-real solution based on finding solutions for the polynomials of order $t$ as truncations of the power series of the $\cos()$ observing the characteristics when $t \to \infty$ . The heuristics suggest that we approach a formal power series whose coefficients blow up without bound (as some not yet estimated function of $t$) when increasing $t$ except the first terms which decreases, but I've no idea of the limit somewhere between zero and $0.5$
So it is surely hopeless to assume a meaningful solution this way.
Here is the attempt to a solution; it might illustrate the occuring problems well.
Let $$f_t(x) = \sum_{k=0}^{t-1} c_k x^k = \sum_{k=0}^{t-1} { ( (î x)^k + (-î x)^k) \over 2 \cdot k!} $$ the polynomial of degree $t-1$ of the $t$ leading terms of the power series for $\cos(x)$. Then we seek the polynomial $$ g_t(x) = \sum_{k=0}^{t-1} a_k { x^k } $$ such that $$g_t(g_t(x)) = f_t(x) + O(x^{t})$$
This process is interesting because in the case of the half-iterate of the $\exp()$ function it seems very likely that this process approximates well the famous Kneser's real-to-real solution (which was mentioned here in MO too). The machinery in Kneser's solution is highly intransparent and Kneser himself did not give an explicite way how to find the power series, however participants in the tetrationforum developed such explicite solutions (or at least asymptotic approximations) which give explicite power series to arbitrary many terms and arbitrary precision.
I found an iterative method to approximate the coefficents in $g_t(x)$ for each $t$ to arbitrary accuracy. The basic principle is the Newton-rootfinding algorithm applied to the (truncated) Carlemanmatrix $F_t$ assigned to $f_t(x)$ finding the (truncated) Carlemanmatrix $G_t$ and from this the assigned function $g_t(x)$ which gives indeed $G_t^2 = \hat F_t$ (where $\hat F_t$ is no more Carleman) . The key ingredient is here, that the Newton-iteration has a restriction to make sure, $G_t$ becomes a true (truncated) Carlemanmatrix - so we might introduce the name "restricted Newton squarerroot finding algorithm on Carlemanmatrices" (I have explained this a bit more in a recent MSE-answer on the half-iterate of the $\exp()$ where also the Kneser-solution was posted.)
The results are the following polynomials $g_t(x)$ which produce perfectly $g_t(g_t(x)) = f_t(x) + O(x^{t})$ that means they reproduce perfectly the $t$-leading coefficients of the $\cos()$-function. Here are the coefficients for the odd $t$ from $t=3$ to $t=21$ (columnwise):
x t=5 t=7 t=9 t=11 t=13 t=15 t=17 t=19 t=21
--+-----------------------------------------------------------------------------------------------------------------
0 0.71233691 0.69301041 0.67288261 0.65596547 0.64204889 0.63051446 0.62082937 0.61258889 0.60549199
1 1.6102585 2.7287951 4.0085148 5.3987263 6.8710340 8.4067966 9.9930135 11.620254 13.281463
2 -3.5729667 -10.358626 -21.756970 -38.394902 -60.724230 -89.082264 -123.72870 -164.86856 -212.66714
3 3.6948540 21.217052 67.801784 162.22728 325.77604 581.54895 954.00757 1468.6586 2151.8298
4 -1.4832464 -24.599214 -132.45279 -450.65133 -1181.2198 -2612.5304 -5126.5486 -9204.2529 -15429.785
5 . 15.320488 166.15792 860.15843 3049.8356 8544.1701 20359.207 43097.180 83346.945
6 . -4.0249864 -130.81135 -1142.7490 -5750.4976 -20980.263 -61813.612 -156251.66 -351914.80
7 . . 59.142056 1043.7623 7979.3068 39297.326 146369.32 448528.28 1189365.8
8 . . -11.778174 -627.54139 -8088.3168 -56426.961 -273197.44 -1033454.9 -3267909.5
9 . . . 224.39055 5842.4420 61818.541 403310.68 1925686.7 7371868.9
10 . . . -36.268480 -2855.3854 -50872.280 -469386.59 -2908877.0 -13728482.
11 . . . . 848.08908 30498.821 426187.16 3553907.6 21143786.
12 . . . . -115.82787 -12594.714 -295982.27 -3485952.6 -26885414.
13 . . . . . 3207.6064 152012.62 2708168.3 28072241.
14 . . . . . -380.22736 -54457.346 -1629854.6 -23835563.
15 . . . . . . 12160.303 733291.62 16205217.
16 . . . . . . -1275.3747 -232281.71 -8615795.2
17 . . . . . . . 46234.919 3452645.4
18 . . . . . . . -4352.9127 -981161.26
19 . . . . . . . . 176317.73
20 . . . . . . . . -15070.867
The coefficients in the polynomials show a clear growth with the degree $t$ and also suggest, that the "naive" extrapolation of the final series would in the leading terms look roughly like a geometric series with some function of $t$ as quotient.
Of course such an extrapolated series is divergent for all $|x|>0$ , and I assume, that a Kneser-like solution is accordingly impossible.
That polynomials $g_t(x)$ are also actually not very useful; while they reproduce the leading terms of the $\cos()$ well, the remaining data in $g_t(g_t(x))$ is much garbage. Here is an example for $t=5$ where $g_5(g_5(x)) = f_5(x) + O(x^5)$ but the $O(x^5)$ -part is really large (and grows with higher polynomials degrees $t$):
g_5(g_5(x)) = f_5(x) + O(x^5) =
1.0000000
-2.2E-44 * x // nonzero because of stopping the Newton-iterations
-0.50000000 * x^2
+2.2E-43 * x^3 // nonzero because of stopping the Newton-iterations
+0.0416667 * x^4
------------------------
+46.474309 * x^5
-292.63771 * x^6
+946.43908 * x^7
-2017.5754 * x^8
+3098.6620 * x^9
-3562.7024 * x^10
+3107.2484 * x^11
-2045.3815 * x^12
+992.01697 * x^13
-336.46574 * x^14
+71.533667 * x^15
-7.1790424 * x^16