An algorithm to compute $\omega$-primality in a numerical monoid