This is discussed a paper by Taussky and Zassenhaus: Symmetric Similarity Transform.
I spend some time understanding this paper, and wrote a note how I understood it in my words. To understand the paper fully, it requires some knowledge on Module theory such as 'Structure Theorem for Finitely Generated Modules over PID'. The following is a starting point:
Case 1: Minimal Polynomial Coincides with Characteristic Polynomial
Let $A$ be a $n\times n$ matrix over a field $K$. Suppose also that the minimal polynomial and characteristic polynomial of $A$ coincide. Then any invertible matrix $X$ satisfying $XA=A^TX$ is symmetric.
Proof)
Consider the following system ($\Sigma_A$) of matrix equations.
$$(1) \ \ XA=A^TX, \ \ (2) \ \ X=X^T.$$
Note that the below system is equivalent to ($\Sigma_A$).
$$(1') \ \ XA=A^TX^T,\ \ (2) \ \ X=X^T.$$
The linear transform $X\mapsto (XA-A^TX^T, \ \ X-X^T)$ has rank at most $n^2-n$ since both components are skew-symmetric.
Thus, the solution space of the system ($\Sigma_A$) has dimension at least $n$.
Now, fix a non-singular transform $X_0$ such that $X_0A=A^TX_0$. Then
$$XA=A^TX \textrm{ if and only if } X_0^{-1}XA=AX_0^{-1}X.$$
This yields an isomorphism $X\mapsto X_0^{-1}X$ between $\{X\mid XA=A^TX\}$ and $C_A=\{X'\mid X'A=AX'\}=\{f(A)|f\in K[t]\}$. Since $\textrm{dim}_{K}C_A = n$, the solution space for (1) has dimension $n$. Then the solution space for ($\Sigma_A$) has dimension at most $n$ because it involves more equations. Therefore, the dimension must be exactly $n$.
Hence, every matrix $X$ satisfying (1) must also satisfy (2).
The following is not discussed in detail in the paper, so I will include it.
Case 2: General Case
Let $A$ be a $n\times n$ matrix over a field $K$. Then there exists an invertible symmetric matrix $X$ such that $XA=A^T X$.
Proof)
Let $X_0$ be an invertible matrix over $K$ such that $A=X_0^{-1} J X_0$ with $J$ is consisted of blocks on diagonal, each block is a companion matrix of $p^s$ for some irreducible polynomial $p$ and $s\geq 1$. Say, $J=\textrm{diag}\{C_1,\ldots, C_r\}$. Now, we may have distinct blocks in $J$ corresponding to the same irreducible polynomial. Consider
$$
XA=A^T X \Longleftrightarrow X_0^{-T}XX_0^{-1}J=J^T X_0^{-T}XX_0^{-1},
$$
where I used the notation $X_0^{-T} = (X_0^{-1})^T$.
By Case 1, we can find a symmetric invertible matrix corresponding to each companion matrix in $J$. Put $ X_0^{-T}XX_0^{-1}=\textrm{diag}\{Y_1,\ldots ,Y_r\}$ such that $Y_i C_i = C_i^T Y_i$ and $Y_i$ is symmetric invertible and of the same size with $C_i$ for each $1\leq i\leq r$.
Then we can put a symmetric invertible matrix
$$
X=X_0^T \textrm{diag}\{Y_1,\ldots ,Y_r\} X_0,
$$
which satisfies $XA=A^T X$.
