How to solve this ODE: $y^{(y(x))}(x)=f(x)$? $$\large{\text{Introduction:}}$$
This question will be partly inspired from:

Evaluation of $$y’=x^y,y’=y^x$$

but what if we made the order of an differential equation equal to the function? Imagine that we had the following linear ordinary differential equation using nth derivative notation. In other words, the $y(x)\,th$ derivative of $y(x)$ is $f(x)$. I would have made $f(x)\to f(x,y)$, but that is too hard to solve. Note that $f(x)$ is any continuous infinitely differentiable function:
$$\frac{d^{y(x)}}{dx^{y(x)}}y(x)=\text D^{y(x)}_x y(x)=f(x)\implies y^{(y)}=f(x)\implies y^{(y)}=f$$
$$\large{\text{Specific Values:}}$$
Here are some examples of points:
$$y^{(y(0))}(0)=f(0), y^{\left(y \left(\frac12\right)\right)} \left(\frac12\right)=f\left(\frac12\right) $$
we could also imagine an inverse relation for $y(x)$ called $y^{-1}(x)$:
$$0=f\left(y^{-1}(0)\right), y^{\left(\frac12\right)} \left(\frac12\right) = f\left(y^{-1}\left(\frac12\right) \right),y’(1)=f\left(y^{-1}(1)\right)$$
$$\large{\text{Problem Statement:}}$$
this presents the problem that if $y(x)\not\in\Bbb Z$, then we need a fractional derivative, but definitions for such derivatives are confusing, so we can use the inverse operation of the fractional derivative which is the fractional integral using the corresponding notation in the bolded link:
$$y^{(y)}=f(x)\implies \text I_x^{y(x)}\text D_x^{y(x)}y(x)=\text I_x^{y(x)}f(x)\implies y+c_1=\text I^y_x f(x)$$
$$\large{\text{Integral Method:}}$$
Now let’s use Cauchy’s Formula for Repeated Integration as the definition for the fractional integration as seen in the bolded link:
$$\,_a\text I^n_x g(x)=\frac 1{Γ(n)}\int_a^xg(t)(x-t)^{n-1}dt\implies y+c_1=\,_a \text I^{y(x)}_x f(x)=\frac{1}{Γ(y(x))}\int_a^x f(t)(x-t)^{y(x)-1} dt\implies \boxed{y!+c_0 Γ(y)=\int_{c_1}^x f(t)(x-t)^{y-1} dt}$$
$$\large{\text{Special Case Solution Using Integral method:}}$$
Let’s now solve the $f(x)=1$ case using the conjectured formula:
$$y^{(y(x))}(x)=1\implies (y(x))!+c_0 Γ(y(x))=\int_{c_1}^x (x-t)^{y(x)-1} dt\implies y!+c_0 Γ(y)=\frac{(x-c_1)^y}y\implies \boxed{y!(y+c_0)-(x-c_1)^y=0}\implies x=\sqrt[y]{y!(y+c_0)}+c_1$$
Here is a complete interactive graph for this conjectured solution.
Another problem is if the fractional derivative has the same definition using $2$ different operator definitions, like this one.
$$\large{\text{Special Case Using Induction:}}$$
Another way for $f(x)=1$ is the following setting all constants of integration to be $0$ via quick induction:
$$y^{(y(x))}(x)=1\mathop\implies^{n=y(x)} y^{(n)}(x)=1\implies y^{(n-1)}=x,y^{(n-2)}=\frac {x^2}{1\cdot2},y^{(n-3)}=\frac{x^3}{3\cdot2\cdot1}\implies y^{(n-k)}=\frac{x^k}{k!}$$
Now let’s set $n=k$:
$$y^{(n-k)}(x)=y^{(n-n)}(x)=y^{(0)}(x)=y(x)=\frac{x^{y(x)}}{(y(x))!}$$
Let’s try to solve for $x$:
$$yx^{-y}=y!$$
but we cannot use the W-Lambert function here. We can also write the conjectured solution as:
$$y^{(y)}=1\mathop\implies^? y\ln(x)-\ln(y)-\ln(y!)=0$$
Using an Inversion theorem will be cumbersome.
Amazingly, we get that the conjectured solution to the differential equation:
$$y^{(y(x))}(x)=1\implies x=\sqrt[y]{yy!}$$
which can be thought of the inverse function of $$y=\sqrt[x]{xx!}$$
graphically shown here
$$\large{\text{Conclusion:}}$$
Is this a correct way to solve $y^{(y)}=f(x)$ and if not, then how? Please correct me and give me feedback!
$$\large{\text{Addendum:}}$$
Also see

Solving $y^{(x)}(x)=ax+b$ in closed form

for a twin question where I was able to find a solution using the following method and assuming all integration constants and $f(x)=1$ for $y(x)\in\Bbb N$:
$$y^{(y(x))}(x)=1=y^{(y)}=1\implies y^{(y-1)}=x+1\implies y^{(x-2)}=\frac{x^2}{2\cdot 1}+x+1$$
Therefore:
$$y^{(y-k)}=\sum_{n=0}^k\frac{x^n}{n!} \implies y^{(y-y)}=y^{(0)}=y(x)=\sum_{n=0}^{y(x)}\frac{x^n}{n!}$$
Here is an interactive graph of the result.
The apparent closed form of the sum may not work also as seen in the question:
$$e^zQ(n,z)=\sum_{m=0}^{n-1}\frac{x^m}{m!},n\in\Bbb N\implies y=\sum_{n=0}^{y} \frac{x^n}{n!}\mathop=^? e^x Q(y+1,x) $$
 A: Not a complete answer, but a partial answer.
I've asked myself the same question before, but never found a complete solution, so I'd like to share my work on it in the hope that it might help someone. Since nobody has answered yet, I will also share the unfinished one. I hope that's ok.
My Work (the best of it)
Let's assume that $f$ has a Maclaurin Series of the general form $f\left( x \right) = \sum_{k = 0}^{\infty}\left( f_{k} \cdot \frac{x^{k}}{k!} \right)$ (Assumption $1.$) and that $f$ can be written as a Generalized Hypergeometric Function (I will abbreviate "Generalized Hypergeometric Function" in the following with the abbreviation "GHF".) of the form $f\left( x \right) = \operatorname{_{p}F_{q}}\left( \begin{matrix} A_{1},\, A_{2},\, \dots,\, A_{p}\\ B_{1},\, B_{2},\, \dots,\, B_{q} \end{matrix};\, x \right)$ (Assumption $2.$), then we could try the Frobenius Method. But before we do that, I would like to mention the nth derivative formula or nth integral formula for monomials that I used:
$$
\begin{align*}
\operatorname{D}_{x}^{n}~ \left( x^{m} \right) &= \frac{m!}{\left( m - n \right)!} \cdot x^{m - n}\\
\operatorname{I}_{x}^{n}~ \left( x^{m} \right) &= \frac{m!}{\left( m + n \right)!} \cdot x^{m + n}\\
\end{align*}
$$
Rewrite it a bit:
$$
\begin{align*}
\operatorname{D}_{x}^{y\left( x \right)} \left( y\left( x \right) \right) &= f\left( x \right)\\
y\left( x \right) &= \operatorname{I}_{x}^{n}~ \left( f\left( x \right) \right)\\
\end{align*}
$$
Substitute $f\left( x \right) = \sum_{k = 0}^{\infty}\left( f_{k} \cdot \frac{x^{k}}{k!} \right)$:
$$
\begin{align*}
y\left( x \right) &= \operatorname{I}_{x}^{n}~ \left( f\left( x \right) \right)\\
y\left( x \right) &= \operatorname{I}_{x}^{n}~ \left( \sum_{k = 0}^{\infty}\left( f_{k} \cdot \frac{x^{k}}{k!} \right) \right)\\
y\left( x \right) &= \sum_{k = 0}^{\infty}\left( f_{k} \cdot \frac{\frac{k!}{\left( k + y\left( x \right) \right)!} \cdot x^{k + y\left( x \right)}}{k!} \right)\\
y\left( x \right) &= \sum_{k = 0}^{\infty}\left( f_{k} \cdot \frac{x^{k + y\left( x \right)}}{\left( k + y\left( x \right) \right)!} \right)\\
y\left( x \right) &= \sum_{k = 0}^{\infty}\left( \frac{f_{k} \cdot x^{y\left( x \right)}}{\left( k + y\left( x \right) \right)!} \cdot x^{k} \right)\\
\end{align*}
$$
We can rewrite this to the generalized series expansion $\sum_{k = 0}^{\infty}\left( c_{k} \cdot \frac{x^{k}}{k!} \right)$ (General Hypergeometric Series) of the GHF:
$$
\begin{align*}
y\left( x \right) &= \sum_{k = 0}^{\infty}\left( \frac{f_{k} \cdot x^{y\left( x \right)}}{\left( k + y\left( x \right) \right)!} \cdot x^{k} \right)\\
y\left( x \right) &= \sum_{k = 0}^{\infty}\left( \frac{f_{k} \cdot x^{y\left( x \right)} \cdot k!}{\left( k + y\left( x \right) \right)!} \cdot \frac{x^{k}}{k!} \right)\\
y\left( x \right) &= \sum_{k = 0}^{\infty}\left( c_{k} \cdot \frac{x^{k}}{k!} \right)\\
\end{align*}
$$
If we try to simplify this series to a GHF, we can get:
$$
\begin{align*}
\frac{c_{k + 1}}{c_{k}} &= \frac{\frac{f_{k + 1} \cdot x^{y\left( x \right)} \cdot \left( k + 1 \right)!}{\left( k + 1 + y\left( x \right) \right)!}}{\frac{f_{k} \cdot x^{y\left( x \right)} \cdot k!}{\left( k + y\left( x \right) \right)!}}\\
\frac{c_{k + 1}}{c_{k}} &= \frac{f_{k + 1} \cdot x^{y\left( x \right)} \cdot \left( k + 1 \right)! \cdot \left( k + y\left( x \right) \right)!}{f_{k} \cdot x^{y\left( x \right)} \cdot k! \cdot \left( k + 1 + y\left( x \right) \right)!}\\
\frac{c_{k + 1}}{c_{k}} &= \frac{f_{k + 1} \cdot \left( k + 1 \right)}{f_{k} \cdot \left( k + 1 + y\left( x \right) \right)}\\
\frac{c_{k + 1}}{c_{k}} &= \frac{f_{k + 1}}{f_{k}} \cdot \frac{k + 1}{k + 1 + y\left( x \right)}\\
\end{align*}
$$
Via the Assumption $2.$ that $f$ can be represented as a GHF of the form $f\left( x \right) = \operatorname{_{p}F_{q}}\left( \begin{matrix} A_{1},\, A_{2},\, \dots,\, A_{p}\\ B_{1},\, B_{2},\, \dots,\, B_{q} \end{matrix};\, x \right)$ it follows:
$$
\begin{align*}
\frac{f_{k + 1}}{f_{k}} &= \frac{\left( k + A_{1} \right) \cdots \left( k + A_{p} \right)}{\left( k + B_{1} \right) \cdots \left( k + B_{q} \right) \cdot \left( k + 1 \right)}\\
f\left( x \right) &= \operatorname{_{p}F_{q}}\left( \begin{matrix}
A_{1},\, A_{2},\, \dots,\, A_{p}\\ B_{1},\, B_{2},\, \dots,\, B_{q}
\end{matrix};\, x \right)\\
\end{align*}
$$
$$
\begin{align*}
\frac{c_{k + 1}}{c_{k}} &= \frac{f_{k + 1}}{f_{k}} \cdot \frac{k + 1}{k + 1 + y\left( x \right)}\\
\frac{c_{k + 1}}{c_{k}} &= \frac{\left( k + A_{1} \right) \cdots \left( k + A_{p} \right)}{\left( k + B_{1} \right) \cdots \left( k + B_{q} \right) \cdot \left( k + 1 \right)} \cdot \frac{k + 1}{k + 1 + y\left( x \right)}\\
\frac{c_{k + 1}}{c_{k}} &= \frac{\left( k + A_{1} \right) \cdots \left( k + A_{p} \right) \cdot \left( k + 1 \right)}{\left( k + B_{1} \right) \cdots \left( k + B_{q} \right) \cdot \left( k + 1 + y\left( x \right) \right) \cdot \left( k + 1 \right)}\\
\end{align*}
$$
We can combine all of this into a GHF:
$$
\begin{align*}
y\left( x \right) &= \operatorname{_{p + 1}F_{q + 1}}\left( \begin{matrix}
1,\, A_{1},\, A_{2},\, \dots,\, A_{p}\\ y\left( x \right) + 1,\, B_{1},\, B_{2},\, \dots,\, B_{q}
\end{matrix};\, x \right)\\
\end{align*}
$$
Leonard Euler found the following Integral Transformation for this form of the GHF:
$$
\begin{align*}
\operatorname{_{A + 1}F_{B + 1}}\left[ \begin{matrix} c,\, a_{1},\, \dots,\, a_{A}\\ d,\, b_{1},\, \dots,\, b_{B}\\ \end{matrix};\, z \right] &= \frac{\Gamma\left( d \right)}{\Gamma\left( c \right) \cdot \Gamma\left( d - c \right)} \cdot \int_{0}^{1} t^{c - 1} \cdot \left( 1 - t \right)^{d - c - 1}\ \cdot \operatorname{_{A}F_{B}}\left[ \begin{matrix} a_{1},\, \dots,\, a_{A}\\ b_{1},\, \dots,\, b_{B}\\ \end{matrix};\, t \cdot z \right]\, \operatorname{d}t\\
\end{align*}
$$
With this we get:
$$
\begin{align*}
y\left( x \right) &= \frac{\Gamma\left( y\left( x \right) + 1 \right)}{\Gamma\left( 1 \right) \cdot \Gamma\left( y\left( x \right) + 1 - 1 \right)} \cdot \int_{0}^{1} t^{1 - 1} \cdot \left( 1 - t \right)^{y\left( x \right) + 1 - 1 - 1}\ \cdot \operatorname{_{p}F_{q}}\left[ \begin{matrix} A_{1},\, \dots,\, A_{p}\\ B_{1},\, \dots,\, B_{q}\\ \end{matrix};\, t \cdot x \right]\, \operatorname{d}t\\
y\left( x \right) &= \frac{\Gamma\left( y\left( x \right) + 1 \right)}{\Gamma\left( 1 \right) \cdot \Gamma\left( y\left( x \right) \right)} \cdot \int_{0}^{1} \left( 1 - t \right)^{y\left( x \right) - 1}\ \cdot \operatorname{_{p}F_{q}}\left[ \begin{matrix} A_{1},\, \dots,\, A_{p}\\ B_{1},\, \dots,\, B_{q}\\ \end{matrix};\, t \cdot x \right]\, \operatorname{d}t\\
y\left( x \right) &= \frac{y\left( x \right)}{\Gamma\left( 1 \right)} \cdot \int_{0}^{1} \left( 1 - t \right)^{y\left( x \right) - 1}\ \cdot \operatorname{_{p}F_{q}}\left[ \begin{matrix} A_{1},\, \dots,\, A_{p}\\ B_{1},\, \dots,\, B_{q}\\ \end{matrix};\, t \cdot x \right]\, \operatorname{d}t\\
\Gamma\left( 1 \right) &= \int_{0}^{1} \left( 1 - t \right)^{y\left( x \right) - 1}\ \cdot \operatorname{_{p}F_{q}}\left[ \begin{matrix} A_{1},\, \dots,\, A_{p}\\ B_{1},\, \dots,\, B_{q}\\ \end{matrix};\, t \cdot x \right]\, \operatorname{d}t\\
1 &= \int_{0}^{1} \left( 1 - t \right)^{y\left( x \right) - 1}\ \cdot \operatorname{_{p}F_{q}}\left[ \begin{matrix} A_{1},\, \dots,\, A_{p}\\ B_{1},\, \dots,\, B_{q}\\ \end{matrix};\, t \cdot x \right]\, \operatorname{d}t\\
\end{align*}
$$
The best result I got was this:
$$\boxed{1 = \int_{0}^{1} \left( 1 - t \right)^{y\left( x \right) - 1}\ \cdot \operatorname{_{p}F_{q}}\left[ \begin{matrix} A_{1},\, \dots,\, A_{p}\\ B_{1},\, \dots,\, B_{q}\\ \end{matrix};\, t \cdot x \right]\, \operatorname{d}t}$$
I don't know if everything I've done is correct, but I hope it can help or at least give an idea for a solution.
