Maxwellian electrodynamics fails when quantum mechanical phenomena are involved, in the same way that Newtonian mechanics needs to be replaced in that regime by quantum mechanics. Maxwell's equations don't really "fail", as there is still an equivalent version in QM, it's just the mechanics itself that changes.
Let me elaborate on that one for a bit. In Newtonian mechanics, you had a time-dependent position and momentum, $x(t)$ and $p(t)$ for your particle. In quantum mechanics, the dynamical state is transferred to the quantum state $\psi$, whose closest classical analogue is a probability density in phase space in Liouvillian mechanics. There are two different "pictures" in quantum mechanics, which are exactly equivalent.
In the Schrödinger picture, the dynamical evolution is encoded in the quantum state $\psi$, which evolves in time according to the Schrödinger equation. The position and momentum are replaced by static operators $\hat x$ and $\hat p$ which act on $\psi$; this action can be used to find the expected value, and other statistics, of any measurement of position or momentum.
In the Heisenberg picture, the quantum state is fixed, and it is the operators of all the dynamical variables, including position and momentum, that evolve in time, via the Heisenberg equation.
In the simplest version of quantum electrodynamics, and in particular when no relativistic phenomena are involved, Maxwell's equations continue to hold: they are exactly the Heisenberg equations on the electric and magnetic fields, which are now operators on the system's state $\psi$. Thus, you're formally still "using" the Maxwell equations, but this is rather misleading as the mechanics around it is completely different. (Also, you tend to work on the Schrödinger picture, but that's beside the point.)
This regime is used to describe experiments that require the field itself to be quantized, such as Hong-Ou-Mandel interferometry or experiments where the field is measurably entangled with matter. There is also a large gray area of experiments which are usefully described with this formalism but do not actually require a quantized EM field, such as the examples mentioned by Anna. (Thus, for example, black-body radiation can be explained equally well with discrete energy levels on the emitters rather than the radiation.)
This regime was, until recently, pretty much confined to optical physics, so it wasn't really something an electrical engineer would need to worry about. That has begun to change with the introduction of circuit QED, which is the study of superconducting circuits which exhibit quantum behaviour. This is an exciting new research field and it's one of our best bets for building a quantum computer (or, depending on who you ask, the model used by the one quantum computer that's already built. ish.), so it's something to look at if you're looking at career options ;).
The really crazy stuff comes in when you push electrodynamics into regimes which are both quantum and relativistic (where "relativistic" means that the frequency $\nu$ of the EM radiation is bigger than $c^2/h$ times the mass of all relevant material particles). Here quantum mechanics also changes, and becomes what's known as quantum field theory, and this introduces a number of different phenomena. In particular, the number of particles may change over time, so you can put a photon in a box and come back to find an electron and a positron (which wouldn't happen in classical EM).
Again, here the problem is not EM itself, but rather the mechanics around it. QFT is built around a concept called the action, which completely determines the dynamics. You can also build classical mechanics around the action, and the action for quantum electrodynamics is formally identical to that of classical electrodynamics.
This regime includes pair creation and annihilation phenomena, and also things like photon-photon scattering, which do seem at odds with classical EM. For example, you can produce two gamma-ray beams and make them intersect, and they will scatter off each other slightly. This is inconsistent with the superposition principle of classical EM, as it breaks linearity, so you could say that the Maxwell equations have failed - but, as I pointed out, it's a bit more subtle than that.
Best Answer
$M$ is a measure of the strength of the dipole field. In SI units, $M= -8\times 10^{15} \ T m^3$, where $T$ stands for "teslas" and $m^3$ means "meters cubed".
If you only want to know what the field looks like, then you could just set $M=-1$, so you have that
$$ B_r = -2\frac{\cos(\theta)}{r^3}$$ $$ B_\theta = -\frac{\sin(\theta)}{r^3}$$ $$ B_\phi = 0$$
As far as actually plotting this goes, you need to be careful. The easiest thing to do would be to convert everything to Cartesian coordinates and then use some kind of vector plotting software, such as the one available through WolframAlpha.
The NASA document gives the cartesian form of the equations on the next page. Keeping our convention that $M=-1$, we have that
$$B_x = -3\frac{xz}{r^5}$$ $$B_y = -3\frac{yz}{r^5}$$ $$B_z = \frac{r^2-3z^2}{r^5}$$
where $r=(x^2+y^2+z^2)^{\frac{1}{2}}$. This is a 3D vector field, which is typically difficult to plot - instead, we can look at a "slice" by simply setting $y=0$. We then take the $z$ axis to be vertical and the $x$ axis to be horizontal.
The components of the vector field are then
$$ B_x = \frac{-3xz}{(x^2+z^2)^{\frac{5}{2}}}$$ $$ B_z = \frac{x^2-2z^2}{(x^2+z^2)^\frac{5}{2}}$$
Plotting this is not as easy at is appears, though, because the quantity $x^2+z^2$ goes to zero at the coordinate origin, which will cause most plotting programs to behave badly.
Additionally, even if you tell the program not to plot any vectors at the origin, the dipole field increases dramatically for smaller $r$. Your plotting program will probably scale the vectors according to the longest ones in the plot, which means you'll end up with a plot that looks like this:
This is not terribly illuminating. What you'll want to do is restrict your attention to $r$ greater than some lower limit. Playing around with it a bit, I think it looks nice if you make the restriction $x^2+z^2>1$ and plot it from $-2$ to $2$ in $x$ and $z$, but you can make that decision yourself.
A quick example of Mathematica syntax:
So a basic plot looks like this:
If I wanted to restrict my vector plot to the region outside the unit circle, I would replace xComponent with
which is equal to "xComponent" if $x^2+y^2>1$, and is equal to $0$ otherwise. That syntax would look like this: