Here the problem is that you are using too many notations and you are changing them constantly. 

1. The master equation in electrostatics : 
$$ \vec{D}= \epsilon_0 \vec{E} + \vec{P}.$$ 

(i) $ \vec{D} $ is known as the **electric displacement vector**. This represents the electric field (But it's really not an true electric field. Check its dimension.) in a system due to the **free charges**. For example, the a charged conductor or ions embedded in a dielectric material are the free charges of this system. We can control these free charges and consequently we have full control over $ \vec{D}$.

(ii) $ \vec{P} $ is the **polarisation vector**. It is defined as the electric dipole moment per unit volume of a system. Polarised the atoms or atoms having permanent dipole moment in the system create tiny dipoles. These tiny dipole moments constitute this polarisation vector. Also note that these tiny dipoles create the **bound charges** in the system. The value of these bound charges can be obtained from the polarisation vector itself. 

(iii) $ \vec{E}$ is the **total electric field** of a system. Means it’s the field due to both free and bound charges present in the system.

(iv) $\epsilon_0$ is obviously the permittivity of free space. Now for **linear dielectrics** (Dielectrics in which polarisation varies linearly with electric field,  $ \vec{E}$.), the defining equation is,
$$ \vec{P} = \epsilon_0 \chi_e  \vec{E}. $$

Note that in RHS we are putting total field $ \vec{E}$ not $ \vec{D}$. This is the definition. By using this convention everything works out well. This also leads to the equation $ \vec{D} = \epsilon \vec{E}  $. Where $ \epsilon $ is the permittivity of the dielectric material. 

Exactly similar quantities appear in magnetostatics too.

2. The master equation in magnetostatics: 
$$ \vec{H} = \vec{B}/\mu_0 - \vec{M}. $$

(i) Here $ \vec{H} $ is known as the **Auxiliary field**. This represents the magnetic field (not true magnetic field) in a system due to the **free currents**. For example the a constant current carrying wire embedded in a paramagnetic material provides the free current to this system. We can control these free currents and consequently we have full control over $ \vec{H}$.

(ii) $ \vec{M} $ is the **magnetisation vector**. It is defined as the magnetic dipole moment per unit volume of a system. The atoms in the system by act as tiny magnetic dipoles. These tiny dipole moments constitute this magnetisation vector. Also note that these tiny dipoles create the **bound currents** in the system. The value of these bound currents can be obtained from the magnetisation vector itself. 

(iii) $ \vec{B}$ is the **total magnetic field** of the system. Means it’s the field due to both free and bound currents present in the system.

(iv) $\mu_0$ is obviously the permeability of free space. Now for **linear magnetic materials** (materials in which magnetisation varies linearly with auxiliary field,  $ \vec{H}$.), the defining equation is,
$$ \vec{M} = \chi_m \vec{H}. $$

This also leads to the equation $ \vec{H} = \vec{B}/\mu  $. Where $ \mu $ is the permeability of the material.  Now let’s come to your problem.

> A sphere of radius R of a linear magnetic material of permeability μ is
> located in a region of empty space where a uniform magnetic field B0
> exists. a) Knowing that the magnetization that appears on the sphere
> is uniform, calculate M, the dipole moment induced in the sphere and
> the B field at all points in space.

Here, auxiliary field $ \vec{H} $ is given as $B_0/\mu_0$ outside the material. Using the equations mentioned above you need to find magnetisation vector $ \vec{M} $ and total magnetic field $ \vec{B}$. As you have mentioned the magnetic field inside the material **due to the constant magnetisation vector only** is $ - \vec{M}/3$. Therefore in this case the total auxiliary field inside the material is $ \vec{H} = (\vec{B_0}/\mu_0) - (\vec{M}/3)$.

So, by putting these in the magnetostatics master equation and using $ \vec{H} = \vec{B}/\mu $ we get,
$$ \frac{\vec{B_0}}{\mu_0} - \frac{\vec{M}}{3} = \frac{\mu}{\mu_0}(\frac{\vec{B_0}}{\mu_0} - \frac{\vec{M}}{3}) - \vec{M}. $$

Simplify the above equation to get your desired result.