As mentioned in [another answer](https://mathoverflow.net/a/421259/124862) the fact that $A$ and $A^T$ are similar corresponds to a commutative diagram of modules over the polynomial ring $F[\lambda]$:
$$
   \begin{array}{c}
       0 & \to & F[\lambda]\otimes V &
         & \stackrel{\lambda I - A}{\longrightarrow} & 
         F[\lambda]\otimes V & \to & V_{A} & \to & 0\\
         && \downarrow~R &&& Q~\downarrow & & \downarrow P \\
       0 & \to & F[\lambda]\otimes V &
         & \stackrel{\lambda I - A^T}{\longrightarrow} &
         F[\lambda]\otimes V & \to & V_{A^T}  & \to & 0     
   \end{array}
$$
where:
  - $V$ denotes that $n$-dimensional column vector space over $F$
  - For a square matrix $B$ we use $V_B$ to denote $V$ as a *module* over $F[\lambda]$ on which $\lambda$ acts by $B$.
  - The maps $P$, $Q$, $R$ are isomorphisms of $F[\lambda]$ modules.

The matrix $P$ is what we need to determine.

The idea is that once we determine $Q$ and $R$, then $P$ can be easily calculated.

The [Wikipedia page for the Smith Normal Form](https://en.wikipedia.org/wiki/Smith_normal_form) sketches an algorithm from which one can compute invertible matrices $F$ and $G$ over $F[\lambda]$ such that $F(\lambda u - A)G$ is a diagonal matrix over $F[\lambda]$. It follows that $G^T(\lambda u -A^T)F^T$ is the *same* diagonal matrix.

We can now take $R^{-1}=G(F^T)^{-1}$ and $Q=(G^{T})^{-1}F$ to get the identity $Q(\lambda I - A)=(\lambda I - A^{T})R$ as required. As mentioned above, now that we have calculated $Q$ and $R$, we can calculate $P$.

**Update**: Since it may not be immediately obvious how to calculate the *matrix* of $P$ with entries in $F$, here is an explicit calculation which uses the matrix $Q$ over $F[\lambda]$.

Note that if $e_i$ denotes the $i$-th standard basis vector of $V$, then under the horizontal map from $F[\lambda]\otimes V$ to $V_A$ or $V_B$, $1\otimes e_i$ goes to $e_i$.

Under $Q$, the image of $1\otimes e_i$ in $F[\lambda]\otimes V$ is $E_i=\sum_{k=1}^n Q_{k,i} \otimes e_k$, where $Q_{k,i}$ is the $(k,i)$-th entry of $Q$ and is a polynomial in $\lambda$. The image of $E_i$ in $V_B$ is $f_i=\sum_{k=1}^n Q_{k,i}(A)e_k$. We then write $f_i=\sum_{l=1} P_{l,i} e_l$ to get the matrix $P$. (It is a bit curious that it only depends on $Q$!) Note that to actually calculate $P$, we need to calculate the matrix entries of $Q_{k,i}(A)$ for all $i$ and $k$.