Without pulling out my EM notes, I think you have a problem in how you are defining your terms as vectors or magnitudes of vectors in the expression.  The first term should have something like a dot product between the dipole and the position vector r.  I suggest an introductory EM text which will walk you through the necessary differential equations.  Griffiths text is quite clear (Introduction to Electrodynamics (4th Edition): David J. Griffiths).  

regards,

David

if you can get it quickly use Wagness  EM fields chpt 19 . Super clear  

or

Use  vector pot. A(r) 

A(r)= { muo * m * sin(theta)/( 4*pi*r^2)  } phi

no Latex in google sorry

where theta & phi are the two   usual spherical polar angles,  m  is dipole moment

B is of course curl( A)  don't forget to use spherical  polar form for the curl  !!

Fairly sure you ca plot  using polars in matlab if not alas tedious trig to get into cartesians

The original expression "((mu)/r^3)*(3(z^2/r^2)-1)" looks like the z component of the magnetic field (i. e., like Bz) of a dipole in z direction (perhaps up to some unit-system-depending factor).

In SI units, the B field (vector) of a dipole m is (cf. link below):

     B(r) = (µ0/4pi) [ 3 r (m · r) / r5  -  m / r3 ]

with r = [x,y,z] and r = (x2+y2+z2)1/2, and for a dipole in z direction we have m = [0,0,mz] and m · r = mz z.

So, this should result in:

      Bx(r) = (µ0/4pi) [ 3 x (mz z) / r5 ]

      By(r) = (µ0/4pi) [ 3 y (mz z) / r5 ]

      Bz(r) = (µ0/4pi) [ 3 z (mz z) / r5 - mz / r3 ]

where the last line is basically your expression from above.

https://en.wikipedia.org/wiki/Magnetic_dipole

Thank you all for the answers. It works now with the original expression !