On the joint numerical range of a pair of symmetric matrices Proposition 13.4 of Alexander Barvinok's A Course in Convexity shows the existence of the following result:

Let $n\ge 3$. For two $n\times n$ symmetric matrices $A$ and $B$, and a PSD matrix $X$ with $\mbox{trace}(X) = 1$, there exists a unit vector $x$ such that $x^T A x = \mbox{trace}(AX)$ and $x^T B x = \mbox{trace}(BX)$.

This book does not show how to construct such a vector!
In Boyd & Vandenberghe's Convex Optimization, on page 656, there exists a constructive method for a looser version of above (basically, no constraints on the trace of $X$ and consequently, the magnitude of $x$).
I could not so far use their proof for the above stronger version result and basically construct such $x$. Any proof, idea, or help?

Some thoughts
Let $\alpha = \mathrm{Tr}(AX)$ and $\beta = \mathrm{Tr}(BX)$.
We need to find $x$ such that
\begin{align}
x^\mathsf{T} A x &= \alpha, \\
x^\mathsf{T} B x &= \beta, \\
x^\mathsf{T} x &= 1.
\end{align}
We use the Newton-Raphson method to solve the above system of three quadratic equations.
The Newton-Raphson map is
$$x \mapsto x - J^\mathsf{T}(JJ^\mathsf{T})^{-1} \left(
                                               \begin{array}{l}
                                                 x^\mathsf{T} A x - \alpha \\
                                                 x^\mathsf{T} B x - \beta \\
                                                 x^\mathsf{T} x - 1
                                               \end{array}
                                             \right)
$$
where $J$ is the Jacobian matrix given by
$$J = \left(
        \begin{array}{l}
          2x^\mathsf{T}A \\
          2x^\mathsf{T}B \\
          2x^\mathsf{T}
        \end{array}
      \right).
$$
Example. Let $n=5$,
$$A = \left(\begin{array}{rrrrr} -6 & -2 & 0 & -1 & 3\\ -2 & 2 & 2 & -2 & -4\\ 0 & 2 & 2 & -5 & 2\\ -1 & -2 & -5 & 4 & -2\\ 3 & -4 & 2 & -2 & 0 \end{array}\right),$$
$$B = \left(\begin{array}{rrrrr} 2 & -1 & -3 & -2 & 2\\ -1 & -6 & -4 & -1 & 0\\ -3 & -4 & 0 & -5 & -1\\ -2 & -1 & -5 & -6 & -4\\ 2 & 0 & -1 & -4 & 6 \end{array}\right)$$
and
$$X = \frac{1}{127}\left(\begin{array}{rrrrr} 29 & -9 & -20 & 16 & -1\\ -9 & 35 & 9 & -20 & 6\\ -20 & 9 & 26 & -22 & -8\\ 16 & -20 & -22 & 24 & 5\\ -1 & 6 & -8 & 5 & 13 \end{array}\right).$$
Randomly choose an initial $x_0$. After a few iterations, we get $x \approx [0.5829,    0.5326,    0.2344,   -0.4843,    0.2951]^\mathsf{T}$.
