I haven't been able to find either tables of these matrix elements or a formula online, so I decided to provide them here myself.
The following tables give the radial matrix elements of $r^k$ for $k$ from -2 to 4, and for the six lowest-energy states (having $n$ from 1 to 3). The rows and columns are labeled by $nl$. The units are such that the Bohr radius $a_0$ is 1.
$$\begin{array}{c|cccccc}
r^{-2} & 10 & 20 & 21 & 30 & 31 & 32 \\
\hline
10 & 2 & \frac{4 \sqrt{2}}{9} & \frac{2 \sqrt{\frac{2}{3}}}{9} & \frac{7}{12 \sqrt{3}} & \frac{1}{4 \sqrt{6}} & \frac{1}{12 \sqrt{30}} \\
20 & \frac{4 \sqrt{2}}{9} & \frac{1}{4} & 0 & \frac{292 \sqrt{\frac{2}{3}}}{1875} & \frac{16}{625 \sqrt{3}} & -\frac{128}{1875 \sqrt{15}} \\
21 & \frac{2 \sqrt{\frac{2}{3}}}{9} & 0 & \frac{1}{12} & \frac{2 \sqrt{2}}{625} & \frac{64}{1875} & \frac{32}{625 \sqrt{5}} \\
30 & \frac{7}{12 \sqrt{3}} & \frac{292 \sqrt{\frac{2}{3}}}{1875} & \frac{2 \sqrt{2}}{625} & \frac{2}{27} & 0 & 0 \\
31 & \frac{1}{4 \sqrt{6}} & \frac{16}{625 \sqrt{3}} & \frac{64}{1875} & 0 & \frac{2}{81} & 0 \\
32 & \frac{1}{12 \sqrt{30}} & -\frac{128}{1875 \sqrt{15}} & \frac{32}{625 \sqrt{5}} & 0 & 0 & \frac{2}{135} \\
\end{array}$$
$$\begin{array}{c|cccccc}
r^{-1} & 10 & 20 & 21 & 30 & 31 & 32 \\
\hline
10 & 1 & \frac{4 \sqrt{2}}{27} & \frac{8 \sqrt{\frac{2}{3}}}{27} & \frac{\sqrt{3}}{16} & \frac{5}{16 \sqrt{6}} & \frac{\sqrt{\frac{3}{10}}}{16} \\
20 & \frac{4 \sqrt{2}}{27} & \frac{1}{4} & -\frac{1}{4 \sqrt{3}} & \frac{92 \sqrt{6}}{3125} & \frac{64}{3125 \sqrt{3}} & -\frac{448 \sqrt{\frac{3}{5}}}{3125} \\
21 & \frac{8 \sqrt{\frac{2}{3}}}{27} & -\frac{1}{4 \sqrt{3}} & \frac{1}{4} & -\frac{72 \sqrt{2}}{3125} & \frac{192}{3125} & \frac{768}{3125 \sqrt{5}} \\
30 & \frac{\sqrt{3}}{16} & \frac{92 \sqrt{6}}{3125} & -\frac{72 \sqrt{2}}{3125} & \frac{1}{9} & -\frac{1}{9 \sqrt{2}} & \frac{1}{9 \sqrt{10}} \\
31 & \frac{5}{16 \sqrt{6}} & \frac{64}{3125 \sqrt{3}} & \frac{192}{3125} & -\frac{1}{9 \sqrt{2}} & \frac{1}{9} & -\frac{1}{9 \sqrt{5}} \\
32 & \frac{\sqrt{\frac{3}{10}}}{16} & -\frac{448 \sqrt{\frac{3}{5}}}{3125} & \frac{768}{3125 \sqrt{5}} & \frac{1}{9 \sqrt{10}} & -\frac{1}{9 \sqrt{5}} & \frac{1}{9} \\
\end{array}$$
$$\begin{array}{c|cccccc}
r^0 & 10 & 20 & 21 & 30 & 31 & 32 \\
\hline
10 & 1 & 0 & \frac{16 \sqrt{\frac{2}{3}}}{27} & 0 & \frac{3 \sqrt{\frac{3}{2}}}{16} & \frac{3 \sqrt{\frac{3}{10}}}{16} \\
20 & 0 & 1 & -\frac{\sqrt{3}}{2} & 0 & \frac{384 \sqrt{3}}{3125} & -\frac{3072 \sqrt{\frac{3}{5}}}{3125} \\
21 & \frac{16 \sqrt{\frac{2}{3}}}{27} & -\frac{\sqrt{3}}{2} & 1 & -\frac{144 \sqrt{2}}{3125} & 0 & \frac{4608}{3125 \sqrt{5}} \\
30 & 0 & 0 & -\frac{144 \sqrt{2}}{3125} & 1 & -\frac{2 \sqrt{2}}{3} & \sqrt{\frac{2}{5}} \\
31 & \frac{3 \sqrt{\frac{3}{2}}}{16} & \frac{384 \sqrt{3}}{3125} & 0 & -\frac{2 \sqrt{2}}{3} & 1 & -\frac{\sqrt{5}}{3} \\
32 & \frac{3 \sqrt{\frac{3}{10}}}{16} & -\frac{3072 \sqrt{\frac{3}{5}}}{3125} & \frac{4608}{3125 \sqrt{5}} & \sqrt{\frac{2}{5}} & -\frac{\sqrt{5}}{3} & 1 \\
\end{array}$$
$$\begin{array}{c|cccccc}
r & 10 & 20 & 21 & 30 & 31 & 32 \\
\hline
10 & \frac{3}{2} & -\frac{32 \sqrt{2}}{81} & \frac{128 \sqrt{\frac{2}{3}}}{81} & -\frac{9 \sqrt{3}}{64} & \frac{27 \sqrt{\frac{3}{2}}}{64} & \frac{9 \sqrt{\frac{15}{2}}}{64} \\
20 & -\frac{32 \sqrt{2}}{81} & 6 & -3 \sqrt{3} & -\frac{11808 \sqrt{6}}{15625} & \frac{27648 \sqrt{3}}{15625} & -\frac{119808 \sqrt{\frac{3}{5}}}{15625} \\
21 & \frac{128 \sqrt{\frac{2}{3}}}{81} & -3 \sqrt{3} & 5 & \frac{10368 \sqrt{2}}{15625} & -\frac{27648}{15625} & \frac{165888}{15625 \sqrt{5}} \\
30 & -\frac{9 \sqrt{3}}{64} & -\frac{11808 \sqrt{6}}{15625} & \frac{10368 \sqrt{2}}{15625} & \frac{27}{2} & -9 \sqrt{2} & 3 \sqrt{10} \\
31 & \frac{27 \sqrt{\frac{3}{2}}}{64} & \frac{27648 \sqrt{3}}{15625} & -\frac{27648}{15625} & -9 \sqrt{2} & \frac{25}{2} & -\frac{9 \sqrt{5}}{2} \\
32 & \frac{9 \sqrt{\frac{15}{2}}}{64} & -\frac{119808 \sqrt{\frac{3}{5}}}{15625} & \frac{165888}{15625 \sqrt{5}} & 3 \sqrt{10} & -\frac{9 \sqrt{5}}{2} & \frac{21}{2} \\
\end{array}$$
$$\begin{array}{c|cccccc}
r^2 & 10 & 20 & 21 & 30 & 31 & 32 \\
\hline
10 & 3 & -\frac{512 \sqrt{2}}{243} & \frac{1280 \sqrt{\frac{2}{3}}}{243} & -\frac{81 \sqrt{3}}{128} & \frac{135 \sqrt{\frac{3}{2}}}{128} & \frac{81 \sqrt{\frac{15}{2}}}{128} \\
20 & -\frac{512 \sqrt{2}}{243} & 42 & -20 \sqrt{3} & -\frac{953856 \sqrt{6}}{78125} & \frac{1714176 \sqrt{3}}{78125} & -\frac{5308416 \sqrt{\frac{3}{5}}}{78125} \\
21 & \frac{1280 \sqrt{\frac{2}{3}}}{243} & -20 \sqrt{3} & 30 & \frac{1057536 \sqrt{2}}{78125} & -\frac{1990656}{78125} & \frac{6967296}{78125 \sqrt{5}} \\
30 & -\frac{81 \sqrt{3}}{128} & -\frac{953856 \sqrt{6}}{78125} & \frac{1057536 \sqrt{2}}{78125} & 207 & -135 \sqrt{2} & 45 \sqrt{10} \\
31 & \frac{135 \sqrt{\frac{3}{2}}}{128} & \frac{1714176 \sqrt{3}}{78125} & -\frac{1990656}{78125} & -135 \sqrt{2} & 180 & -63 \sqrt{5} \\
32 & \frac{81 \sqrt{\frac{15}{2}}}{128} & -\frac{5308416 \sqrt{\frac{3}{5}}}{78125} & \frac{6967296}{78125 \sqrt{5}} & 45 \sqrt{10} & -63 \sqrt{5} & 126 \\
\end{array}$$
$$\begin{array}{c|cccccc}
r^3 & 10 & 20 & 21 & 30 & 31 & 32 \\
\hline
10 & \frac{15}{2} & -\frac{2560 \sqrt{2}}{243} & \frac{5120 \sqrt{\frac{2}{3}}}{243} & -\frac{1215 \sqrt{3}}{512} & \frac{1215 \sqrt{\frac{3}{2}}}{512} & \frac{1701 \sqrt{\frac{15}{2}}}{512} \\
20 & -\frac{2560 \sqrt{2}}{243} & 330 & -150 \sqrt{3} & -\frac{66313728 \sqrt{6}}{390625} & \frac{105504768 \sqrt{3}}{390625} & -\frac{264757248 \sqrt{\frac{3}{5}}}{390625} \\
21 & \frac{5120 \sqrt{\frac{2}{3}}}{243} & -150 \sqrt{3} & 210 & \frac{76889088 \sqrt{2}}{390625} & -\frac{125411328}{390625} & \frac{334430208}{390625 \sqrt{5}} \\
30 & -\frac{1215 \sqrt{3}}{512} & -\frac{66313728 \sqrt{6}}{390625} & \frac{76889088 \sqrt{2}}{390625} & \frac{6885}{2} & -\frac{8775}{2 \sqrt{2}} & \frac{2835 \sqrt{\frac{5}{2}}}{2} \\
31 & \frac{1215 \sqrt{\frac{3}{2}}}{512} & \frac{105504768 \sqrt{3}}{390625} & -\frac{125411328}{390625} & -\frac{8775}{2 \sqrt{2}} & 2835 & -945 \sqrt{5} \\
32 & \frac{1701 \sqrt{\frac{15}{2}}}{512} & -\frac{264757248 \sqrt{\frac{3}{5}}}{390625} & \frac{334430208}{390625 \sqrt{5}} & \frac{2835 \sqrt{\frac{5}{2}}}{2} & -945 \sqrt{5} & 1701 \\
\end{array}$$
$$\begin{array}{c|cccccc}
r^4 & 10 & 20 & 21 & 30 & 31 & 32 \\
\hline
10 & \frac{45}{2} & -\frac{40960 \sqrt{2}}{729} & \frac{71680 \sqrt{\frac{2}{3}}}{729} & -\frac{3645 \sqrt{3}}{512} & 0 & \frac{5103 \sqrt{\frac{15}{2}}}{256} \\
20 & -\frac{40960 \sqrt{2}}{729} & 2880 & -1260 \sqrt{3} & -\frac{4592443392 \sqrt{6}}{1953125} & \frac{6772211712 \sqrt{3}}{1953125} & -\frac{14714929152 \sqrt{\frac{3}{5}}}{1953125} \\
21 & \frac{71680 \sqrt{\frac{2}{3}}}{729} & -1260 \sqrt{3} & 1680 & \frac{5361334272 \sqrt{2}}{1953125} & -\frac{8026324992}{1953125} & \frac{18059231232}{1953125 \sqrt{5}} \\
30 & -\frac{3645 \sqrt{3}}{512} & -\frac{4592443392 \sqrt{6}}{1953125} & \frac{5361334272 \sqrt{2}}{1953125} & \frac{122715}{2} & -\frac{76545}{\sqrt{2}} & 11907 \sqrt{10} \\
31 & 0 & \frac{6772211712 \sqrt{3}}{1953125} & -\frac{8026324992}{1953125} & -\frac{76545}{\sqrt{2}} & 48195 & -15309 \sqrt{5} \\
32 & \frac{5103 \sqrt{\frac{15}{2}}}{256} & -\frac{14714929152 \sqrt{\frac{3}{5}}}{1953125} & \frac{18059231232}{1953125 \sqrt{5}} & 11907 \sqrt{10} & -15309 \sqrt{5} & 25515 \\
\end{array}$$
To clarify, these radial matrix elements are
$$\langle n' l' | \, r^k \, | n l \rangle \equiv \int_0^\infty r^2 dr \, R_{n'l'}(r) \, r^k \, R_{nl}(r)$$
where
$$R_{nl}(r) = N_{nl} \, e^{-\rho/2} \rho^l L_{n-l+1}^{2l+1}(\rho)$$
is the radial part of the hydrogen wavefunction
$$\psi_{nlm}(r, \theta, \phi) = R_{nl}(r) Y_l^m(\theta, \phi),$$
$\rho$ is related to $r$ by
$$\rho \equiv \frac{2r}{n},$$
and the normalization factor $N_{nl}$ is
$$N_{nl} \equiv \sqrt{\sqrt{\left(\frac{2}{n}\right)^3} \frac{(n-l-1)!}{2n(n+l)!}}.$$
Remember, I am using units where $a_0=1$.
The integral for a particular matrix element can be done in any computer algebra system; it’s just a polynomial times an exponential. The TeX output for the tables above was generated by Mathematica.
A general formula for an arbitrary matrix element can be obtained by using the following integral formula from functions.wolfram.com:
$$I(\alpha, p; m, \lambda, a; n, \beta, b) \equiv
\int_0^\infty t^{\alpha-1} e^{-pt} L_m^\lambda(at) L_n^\beta(bt)\,dt \\
=\frac{\Gamma(\alpha)}{p^\alpha} \frac{(\lambda+1)_m}{m!} \frac{(\beta+1)_n}{n!} \sum_{j=0}^m \frac{(-m)_j(\alpha)_j}{(\lambda+1)_j j!} \left(\frac{a}{p}\right)^j \sum_{k=0}^n \frac{(-n)_k(j+\alpha)_k}{(\beta+1)_k k!} \left(\frac{b}{p}\right)^k$$
where $(x)_n$ is the Pochhammer symbol
$$(x)_n \equiv x(x+1)\dots(x+n-1).$$
The result for the radial matrix element is
$$\langle n' l' | \, r^k \, | n l \rangle = N_{n'l'} N_{nl} \left(\frac{2}{n'}\right)^{l'} \left(\frac{2}{n}\right)^l \\
\times I\left(3+k+l'+l, \frac{1}{n'}+\frac{1}{n}; n'-l'-1, 2l'+1, \frac{2}{n’}; n-l-1, 2l+1,\frac{2}{n}\right).$$
Note that the double-sum $I$ has $(n'-l')(n-l)$ terms, each of which is rational if the $k$ in $r^k$ is. The square roots in the matrix elements come from the normalization factors.
This formula is considerably quicker (at least in Mathematica) than evaluating the individual integrals, especially when $n$ or $n’$ get large.
This formula works for non-integral $k$ as well.
