This is just a speculation so I'm not sure they found the solution this way, but this is 'doable' in my opinion. The idea comes from observing the $x$ that $S^n(x)$ changes its increasing\decreasing state. ($S(x)=4x(1-x)$ in the book)
$$S(1/2)=1\\S^2(\frac{2-\sqrt{2}}{4})=S^2(\frac{2+\sqrt{2}}{4})=1\\S^3(\frac{2-\sqrt{2-\sqrt{2}}}{4})=S^3(\frac{2+\sqrt{2-\sqrt{2}}}{4})=S^3(\frac{2-\sqrt{2+\sqrt{2}}}{4})=S^3(\frac{2+\sqrt{2+\sqrt{2}}}{4})=1\\\vdots$$
Which can be rewritten as
$$S(\sin^2\frac{\pi}{4})=1\\S^2(\sin^2\frac{\pi}{8})=S^2(\sin^2\frac{3\pi}{8})=1\\S^3(\sin^2\frac{\pi}{16})=S^3(\sin^2\frac{3\pi}{16})=S^3(\sin^2\frac{5\pi}{16})=S^3(\sin^2\frac{7\pi}{16})=1\\\vdots$$
Now assume $Pf=f$. From an intuitive viewpoint, the transformation $P$ sends $S^{-1}(x)$ to $x$. So we get the idea to define a function $g : [0,1]\rightarrow [0,1]$ as$$g(x)=f(\sin^2\frac{\pi x}{2})$$
Then the functional equation of $f$ becomes much simpler.
$$g(x)=f(\sin^2\frac{\pi x}{2})=(Pf)(\sin^2\frac{\pi x}{2})=\frac{f(\sin^2\frac{\pi x}{4})+f(\sin^2(\frac{\pi}{2}-\frac{\pi x}{4}))}{4\cos\frac{\pi x}{2}}=\frac{g(\frac{x}{2})+g(1-\frac{x}{2})}{4\cos\frac{\pi x}{2}}$$
It is now easy to see
$$g(\frac{x}{2})\sin\frac{\pi x}{2}+g(1-\frac{x}{2})\sin\pi(1-\frac{x}{2})=4g(x)\sin\frac{\pi x}{2}\cos\frac{\pi x}{2}=2g(x)\sin\pi x$$
Putting $h(x)=g(x)\sin\pi x$, we get $$\frac{1}{2}(h(\frac{x}{2})+h(1-\frac{x}{2}))=h(x)$$
We will now see that the only possible continuous function $h$ is the constant function. Using this equation repetitively, $$h(x)=\frac{1}{2^n}\Big(\sum_{k=1}^{2^{n-1}-1}\big(h(\frac{k}{2^{n-1}}+\frac{x}{2^n})+h(\frac{k}{2^{n-1}}-\frac{x}{2^n})\big)+h(\frac{x}{2^n})+h(1-\frac{x}{2^n})\Big)$$
If x is irrational, each intervals $(\frac{t}{2^n}, \frac{t+1}{2^n})$ contain exactly one of the $\frac{k}{2^{n-1}}\pm\frac{x}{2^n}$ or $\frac{x}{2^n}$ or $1-\frac{x}{2^n}$. (More specifically, $\frac{k}{2^{n-1}}+\frac{x}{2^n}\in(\frac{2k}{2^n}, \frac{2k+1}{2^n})$ and similarly for others.)
So consindering the partition $0, \frac{1}{2^n}, \frac{2}{2^n}, \cdots, 1$ the above is a riemann sum. Because $h$ is a continous function, the sum converges to the Riemann integral which makes $$h(x)=\int_{0}^{1}h(t)dt=c$$ for all irrational $x$. Since the irrational numbers in $[0,1]$ are dense in $[0,1]$, $h(x)=c$ for all $[0,1]$ by the continuity of $h$.
So we finally have $h(x)=c$ and from this we can get $$f(x)=\frac{c}{\sqrt{x(1-x)}}$$
Putting $c=1/\pi$ makes the function $f$ a proabability density function on $[0,1]$.
Update. I've just found out that this was actually a method called the 'change of variables'. It is illustrated Ch6.5 in the OP's book and the idea is to transform $S$ to an exact function $T=g\circ S\circ g^{-1}$ in order to make it more easy to manipulate.(in this case, to the tent function)
According to the book's theorem 6.5.2, if we let $T:(0,1)\rightarrow (0,1)$ a measurable nonsingular transformation, $\phi\in D((a,b))$ a positive density, and $S:(a,b)\rightarrow (a,b)$ defined as $S=g^{-1}\circ S\circ g$ with $$g(x)=\int_{a}^{x}\phi (y)dy$$ then $T$ is exact iff $S$ is statistically stable and $\phi$ is the measure invariant to $S$. So the method only requires to find such $\phi$.
The problem is that this idea was originally due to Ruelle[1977] and Pianigiani[1983] about 30 years later which means that it's not likely for this method that was used by Ulam and Neumann to find $f^*$. So I guess my answer becomes "Well, there is a method for finding such function. But I'm not sure how they found it...".
