According to Dr James Keeler of Cambridge University in one of his lectures at Queensland University, the angle theta can be any angle, but I struggle to see how this is true. There are only a certain number of I values that nuclei take and as illustrated in the figure above, each magnetic quantum number corresponds to a fixed theta. Is my logic flawed?
You are correct that a given value of $m_I$ and $I$ does correspond to one single value of the angle $\theta$. For the most commonly used nuclei in NMR spectroscopy, the "allowed" values are $I = 1/2$ and $m_I = \pm 1/2$. The two states are commonly labelled $\alpha$ and $\beta$; these are the eigenstates of the angular momentum projection operator $\hat{I}_z$, such that
$$\begin{align} \hat{I}_z|\alpha\rangle &= \frac{1}{2}|\alpha\rangle \\ \hat{I}_z|\beta\rangle &= -\frac{1}{2}|\beta\rangle \end{align}$$
(treating $\hbar \equiv 1$). The problem is that you are assuming that the wavefunction (of a spin) can only be in an eigenstate. This is not true. It is perfectly valid, and indeed necessary (see below), to have states that are superpositions of eigenstates. In general, a spin can have any state
$$|\psi\rangle = c_\alpha |\alpha\rangle + c_\beta |\beta\rangle$$
where $|c_\alpha|^2 + |c_\beta|^2 = 1$ (normalisation), and in general the coefficients are time-dependent (see below). This is no longer an eigenstate of $\hat{I}_z$, and consequently is no longer well-described by one value of $m_I$, or for that matter, $\theta$. But that is not intrinsically a problem. The only circumstance in which a wavefunction must be described by an eigenstate is immediately following an observation, where the wavefunction $|\psi\rangle$ collapses to the eigenstate $|i\rangle$ with probability $|c_i|^2$. Otherwise, superposition states are entirely consistent with quantum mechanics.
A more generalised example is to think of the time-independent Schrödinger equation:
$$\hat{H}|\psi\rangle = E|\psi\rangle$$
which is an eigenvalue equation for the Hamiltonian, $\hat{H}$. It is true that only eigenstates of the Hamiltonian are admissible solutions to this equation. But that does not mean that the only physically permissible states for the system are eigenstates of the Hamiltonian. The time-independent SE only governs systems that are not evolving in time, or stationary states. We need something else to describe systems that do evolve with time, which is the time-dependent Schrödinger equation:
$$\hat{H}|\Psi\rangle = i\hbar\frac{\partial |\Psi\rangle}{\partial t}$$
The most general solution to the TDSE is, indeed, a superposition of stationary states (the Internet has lots of more thorough explanations, or read Griffiths' Introduction to Quantum Mechanics):
$$|\Psi\rangle = \sum_i c_i \exp\left(-\frac{iE_it}{\hbar}\right) |\psi_i\rangle $$
where the $c_i$'s are time-independent. Indeed, it is necessary for us to have superposition states, because this allows the coefficients of each eigenstate to evolve in a time-dependent fashion (because of the exponential term, which evolves at different frequencies for different eigenstates). If a system is prepared in an eigenstate of the Hamiltonian, then it can never evolve into another eigenstate (the exponential evolves, but if there is only one eigenstate, it is no more than a simple phase factor).
This might not seem immediately relevant for NMR, but recall that the (free) Hamiltonian in NMR is simply $-\gamma B_0 \hat{I}_z$. So, the eigenstates of $\hat{I}_z$, $|\alpha\rangle$ and $|\beta\rangle$, are also the eigenstates of the Hamiltonian. If we never had superposition eigenstates, then the spins would simply sit there forever and do nothing.*
* As long as the Hamiltonian doesn't change. But if you change the Hamiltonian, say from $\hat{H}_0$ to $\hat{H}_1$, then the eigenstates change, and so your system (which was originally an eigenstate of $\hat{H}_0$ is no longer an eigenstate of $\hat{H}_1$). Indeed, it would be in a superposition of eigenstates of $\hat{H}_1$, and would evolve in time according to the TDSE, so it doesn't invalidate any of the preceding discussion.
The magnetization is composed of the sum of dipole moments associated with individual spins:
$$\vec{M}=\sum_i \vec{\mu}_i$$
and the component along the magnetizing field (assumed to lie along the z axis) is
$$\begin{align}M_z&= \left( \sum_i \vec{\mu}_i\right)_z \\ &= \sum_i \mu_{iz}\end{align}$$
Since the individual dipole moment components along z take on discrete values (are quantized, ie for the homonuclear case $\vec{\mu}_{iz}=\gamma\hbar m_z $) we can write
$$M_z = \gamma\hbar\sum_i m_{i} \tag{1}$$
where the symbols have their usual meanings. Since the sum extends over a large number of spins $i$, the resulting magnetization can for practical purposes be regarded as taking on a continuous range of values, even if the z-component of individual spins is discretely valued.
