Hydrogenic Wavefunctions:
A hydrogenic wavefunction, $\psi_{nlm}(\boldsymbol{r})$, can be written as the product
$$
\psi_{nlm}(r,\theta,\phi) = R_{nl}(r) Y_{lm}(\theta,\phi)
$$
where $R_{nl}$ is a radial wavefunction and $ Y_{lm}$ is a spherical harmonic.
- $R_{nl}$ is a radial wavefunction
- $ Y_{lm}$ is a spherical harmonic
- $n$ is the principal quantum number, $n = 1,2,3\dots$
- $l$ is the angular quantum number, $l = 0=s, 1=p, 2=d,\dots,n-1$
- $m$ is the magnetic quantum number $m = -l, -l+1,\dots,l-1,l$
Probability/Electron Density
The probability density, $\rho$, is equal to the square modulus of the wavefunction:
$$
\begin{align}
\rho_{nlm}(\boldsymbol{r}) &= |\psi_{nlm}(\boldsymbol{r})|^2 \\
\rho_{nlm}(r,\theta,\phi)&= |R_{nl}(r)|^2 \cdot |Y_{lm}(\theta,\phi)|^2
\end{align}
$$
The spherical harmonics $Y_{lm}$ have angular nodes through the nuclei for all values of $l$ and $m$ except $l =m =0$, so we can say that for all $l>0$ orbitals, i.e. all orbitals except s-orbitals, the nuclei will have $0$ probability density.
S-orbitals
Consider s-orbitals, in atomic units and substituting in the solution of the radial wavefunctions in generalized Laguerre polynomials $L^\alpha_k$:
$$
\begin{align}
\rho_{n00}(\boldsymbol{r})&= |R_{n0}(r)|^2 \cdot |Y_{00}(\theta,\phi)|^2\\
&= |R_{n0}(r)|^2 \frac{1}{4\pi} \\
\rho_{n00}(r)&= \frac{Z^3}{\pi n^5} \exp\left(-Zr/n\right) L^1_{n-1}\left(2Zr/n\right) \\
&= A(n) \exp\left(-Zr/n\right) L^1_{n-1}\left(2Zr/n\right)
\end{align}
$$
For all values of $n$, $L^1_{n-1}\left(0\right) = n$, so we know that the nucleus never has zero probability density for s-orbitals, and has probability density
$$
\rho_{n00}(r=0) = \frac{Z^3}{\pi n^4}
$$
Stationary Points
But we can differentiate with respect to $r$ to find stationary points in the radial probability density.
$$
\begin{align}
\frac{d}{dr}\rho_{n00}(r)&= A(n) \frac{d}{dr} \exp\left(-Zr/n\right) L^1_{n-1}\left(2Zr/n\right) \\
&= A(n) \exp\left(-Zr/n\right)\left[ -\frac{Z}{n} L^1_{n-1}\left(2Zr/n\right) + \frac{d}{dr}L^1_{n-1}\left(2Zr/n\right) \right] \\
\text{if}\;n=1: \frac{d}{dr}\rho_{000}(r)&= -\frac{Z^4}{\pi n^6} \exp\left(-Zr/n\right) \\
&> 0 \;\forall\;r
\end{align}
$$
Hence the probability density is constantly decreasing away from the nucleus. So as you previously assumed the maximum for probability density for the 1s orbital is at the nucleus.
$n>1$:
$$
\begin{align}
\text{if}\;n>1 : \frac{d}{dr}\rho_{n00}(r)&= A(n) \exp\left(-Zr/n\right)\left[ -\frac{Z}{n} L^1_{n-1}\left(2Zr/n\right) + \frac{d}{dr}L^1_{n-1}\left(2Zr/n\right) \right] \\
&= A(n) \exp\left(-Zr/n\right)\left[ -\frac{Z}{n} L^1_{n-1}\left(2Zr/n\right) - L^2_{n-2}\left(2Zr/n\right) \right]
\end{align}
$$
Setting the expression above equal to zero and solving for $r$ will give the positions for the radial nodes and maxima inbetween for a given value of $n$. I am not aware of a general method/expression in n.
As you asked about nucleus, we can also consider the point $r=0$:
$$
\begin{align}
\text{if}\;n>1 : \frac{d}{dr}\rho_{n00}(r)&= - \frac{Z^3}{\pi n^5} \left[\frac{Z}{n} L^1_{n-1}\left(0\right) + L^2_{n-2}\left(0\right) \right] \\
&= - \frac{Z^3}{\pi n^5} \left[\frac{Z}{n} n +\frac{n^2-n}{2}\right] \\
&= - \frac{Z^4}{\pi n^5} - \frac{Z^3 \left(n-1\right)}{2\pi n^4} < 0
\end{align}
$$
Conclusion/TL;DR
Hence the nuclei is always a local maximum in the probability density for an s-orbital. You would expect the exponential decay in the radial wavefunction to dominate over the Laguerre polynomial term, so would expect local maxima further out to be smaller that the maxima at the nuclei and hence that the point of highest probability density is at the nucleus.
This makes physical sense, as classically the negatively charged electron wants to be as close as possible to the positively charged nucleus to minimise its potential energy, but the quantum mechanical nature of the electron prevents it from tumbling directly into the nucleus.
