  Given a square matrix $B$ over a principal ideal domain $D$ and an ideal $J$
of $D$, the $J$-ideal of $B$ consists of the polynomials $f\in D[X]$ such that
all entries of $f(B)$ are in $J$. It has been shown that in order to determine
all $J$-ideals of $B$ it suffices to compute a generating set of the
$(p^t)$-ideal of $B$ for finitely many prime powers $p^t$. Moreover, it is
known that $(p^t)$-ideals are generated by polynomials of the form
$p^{t-s}\nu_s$ where $\nu_s$ is a monic polynomial of minimal degree in the
$(p^s)$-ideal of $B$ for some $s\le t$. However, except for the case of
diagonal matrices, it was not known how to determine these polynomials
explicitly. We present an algorithm which allows to compute the polynomials
$\nu_s$ for general square matrices. Exploiting one of McCoy's theorems we
first compute some set of generators of the $(p^s)$-ideal of $B$ which then can
be used to determine $\nu_s$. This algorithmic computation significantly
extends our understanding of the $J$-ideals of $B$.
