You need to look at the idea of Separation of Variables for Partial Differential Equations.
You consider a toy universe comprising oscillators in a box: let's think of a cuboid microwave cavity with electromagnetic fields losslessly confined within perfectly reflecting walls.
The Cartesian components of the electric field all fulfill Helmholtz's equation:
$$\left(\nabla^2 + \frac{\omega^2}{c^2}\right)\psi = 0$$
and they must vanish at all the boundaries. If you work though the standard separation of variables technique, you get a set of allowable modes with standing plane wave variation for the vibration and any general solution is superposition of these allowedd modes. So, the $x$, $y$, and $z$ variations are all of the form $\sin(k_j\,x_j)$.
The electric fields must vanish at the boundaries, this means that the $k_j$ can only be discrete values: $\frac{pi}{L_x},\,\frac{3\,pi}{L_x},\,\frac{5\,pi}{L_x},\,\cdots$ for the $k_x$ values, where we suppose the reflecting walls are at $x=0$ and $x=L_x$. Likewise, the $k_y$ and $k_z$ are restricted to $\frac{(2 n_j + 1)pi}{L_j}$ for $n_j = 0,\,1\,2,\,\cdots$. So we have modes defined by mode numbers: triplets of integers $(n_x,\,n_y,\,n_z)$, and the total variation of each mode is of the form:
$$\psi_{(n_x,\,n_y,\,n_z)}(x,\,y,\,z,\,t) = \sin\left(\frac{(2\,n_x+1)\,x}{L_x}\right)\,\sin\left(\frac{(2\,n_y+1)\,x}{L_y}\right)\,\sin\left(\frac{(2\,n_z+1)\,z}{L_z}\right)\\\,\cos(\omega(n_x,\,n_y,\,n_z)\,t + \delta)$$
But now, substitute this back into the Helmholtz equation and we get:
$$\omega(n_x,\,n_y,\,n_z)^2 = c^2 \left(\frac{(2\,n_x+1)^2}{L_x^2}+\frac{(2\,n_y+1)^2}{L_y^2} + \frac{(2\,n_z+1)^2}{L_z^2}\right)$$
and the above equation can only be fulfilled for discrete $\omega(n_x,\,n_y,\,n_z)$, given that the $n_j$ are integers. Now think of the integer triplets as points in 3D space: the number of modes of frequency $\omega$ or less is the number of these discrete points in the positive $x,\,y\,z$ eighth of the sphere with radius less than $\omega/c$ - in other words, roughly the volume of this eigth, as $\omega$ gets large. Thus the rate of increase of allowed modes with $\omega$ is roughly proportional to the surface area of the spherical sector. In other words, proportional to $\omega^2$. There are roughly $100$ times more modes between frequency $10\,\omega_0$ and $10\,\omega_0 + \Delta$ as there are between frequency $\omega_0$ and $\omega_0 + \Delta$.