The given formula (in that website) is correct.
Assume we have a cylinder with permanent (axial) magnetization $\mathbf{M}$. It will cause a surface current density $\mathbf{K}$:
$$\cases{\bf{K}=\bf M \times \hat n \\ \mathbf M=M_0\hat z}\to \mathbf{K}=M_0\hat \phi$$
So it is like a finite length solenoid. To find the field on it's axis ($z$, point $P$ in the below picture), from Biot–Savart law we'll arrive at:
$$dB=dB_z=\frac{\mu_0K\mathrm{d}z}{2}\frac{R^2}{\xi^2}$$
$$\cases{dz=\frac{Rd\theta}{\sin^2\theta}\\ \frac{R}{z}=\tan \theta \\ \xi=\frac{R}{\sin \theta}}\to dB=-\frac{\mu_0 K}{2}\sin \theta d\theta$$
$$B=\int_{\theta_1}^{\theta_2}dB=\frac{\mu_0 K}{2}(\cos \theta_1-\cos \theta_2)$$
so
$$\boxed {\mathbf{B}=\frac{\mu_0 M_0}{2}(\cos \theta_1-\cos \theta_2)\hat z}$$