It sounds like you're trying to find the shifts in the energy levels caused by $H_V$. To find the energy shift that goes as the first order in your small parameter you compute
$$E_n^{(1)} = \langle \psi_n^{(0)}|H_V|\psi_n^{(0)} \rangle \quad (1)$$
as you noted. In this expression $|\psi_n^{(0)}\rangle$ means "the $n^{\textrm{th}}$ eigenstate of the unperturbed Hamiltonian." That superscript (0) means "to zeroth order in the perturbation" ie. without any perturbation.
So to compute the energy level shifts you have to know the eigenstates of the unperturbed Hamiltonian. You say you already found those to be $|1,\pm 1\rangle$, $|1,0\rangle$ and $|0,0\rangle$. To make contact with the notation in equation (1) let's go ahead and number these
$$|\psi_1^{(0)}\rangle = |0,0 \rangle \quad |\psi_2^{(0)}\rangle = |1,0 \rangle \quad |\psi_3^{(0)}\rangle = |1,-1 \rangle \quad |\psi_4^{(0)}\rangle = |1,+1 \rangle$$
Now you want to compute the energy shifts. Let's just do one of them as an example. Let's compute $E_1^{(1)}$.
$$\begin{eqnarray}
E_1^{(1)} &=& \langle \psi_1^{(0)}|H_V|\psi_1^{(0)}\rangle \\
&=& -\mu B \cdot \langle0,0|\sigma_1 + \sigma_2 | 0,0 \rangle \\
&=& -\mu B \cdot \left( \langle0,0|\sigma_1|0,0\rangle + \langle0,0|\sigma_2|0,0\rangle \right) \end{eqnarray}$$
Now, those two operators $\sigma_1$ and $\sigma_2$ are the spin operators for each individual particle. The state $|0,0\rangle$ is expressed in a basis that is not the single particle basis. $|0,0\rangle$ is a state designed to be simple in the basis of total spin. In order to work out the matrix elements you need for the computation you have to either convert the operators to the total spin basis, or re-express the states in the individual spin basis. I think the latter is easier in this case. Just use
$$|0,0\rangle \equiv \left( \uparrow \downarrow - \downarrow \uparrow \right)/\sqrt{2}.$$
Now it's easy to compute what you need. Here's an example
$$B \cdot \langle 0,0|\sigma_1| 0,0\rangle = \frac{1}{2} \langle\uparrow \downarrow-\downarrow \uparrow|B^x\sigma^x_1 + B^y\sigma^y_1 + B^z\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle.$$
You asked about what the matrices for $\sigma$ should be. Now that we've got everything expressed in the single particle spin basis, it's easy. For example,
$$\sigma^x_1 = \left(\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array} \right).$$
The thing to note is that this only acts on the first particle, so you get things like this
$$\begin{eqnarray}
\langle \uparrow \downarrow | \sigma^x_1 | \downarrow \uparrow \rangle &=& \langle \uparrow | \sigma^x_1 | \downarrow \rangle\cdot\langle \downarrow|\uparrow \rangle \\
&=& \langle \uparrow|\uparrow \rangle \cdot \langle \downarrow|\uparrow \rangle \\
&=& 1 \cdot 0 \\
&=& 0.
\end{eqnarray}$$
At this point I'm sure you can figure out what's going on. I think your problem was one of notation and hopefully seeing this written out helps.
Extension:
I've shown how to compute simple terms like
$$\langle \uparrow \downarrow | \sigma^x_1 | \downarrow \uparrow \rangle$$
but the original expression to compute was more complex. In particular we had
$$\langle\uparrow \downarrow-\downarrow \uparrow|B^x\sigma^x_1 + B^y\sigma^y_1 + B^z\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle.$$
The trick here is that this is all linear so we can break it into little parts like this
$$\langle\uparrow \downarrow-\downarrow \uparrow|B^x\sigma^x_1 + B^y\sigma^y_1 + B^z\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle$$
$$=B_x\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^x_1|\uparrow\downarrow - \downarrow\uparrow\rangle + B_y\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^y_1|\uparrow\downarrow - \downarrow\uparrow\rangle + B_z\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle.$$
Each one of these terms can be broken down further:
$$\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^x_1|\uparrow\downarrow - \downarrow\uparrow\rangle$$
$$\langle \uparrow \downarrow|\sigma^x_1|\uparrow\downarrow\rangle - \langle \uparrow\downarrow | \sigma^x_1 | \downarrow \uparrow\rangle - \langle \downarrow \uparrow | \sigma^x_1 | \uparrow \downarrow \rangle + \langle \downarrow \uparrow | \sigma_1^x | \downarrow \uparrow \rangle$$
From here you can surely see how to to calculate the desired quantity. Note that if you think carefully you can tell that many terms are zero without having to calculate them.
What you do in perturbation theory is you assume the correct eigenvalue equation for a (hitherto) unknown correct wavefuction $|\psi_n\rangle$ and its associated eigenvalue $E_n$: $$ H|\psi_n\rangle = E_n |\psi_n\rangle,$$ where $H$ is the full Hamiltionian.
What you have is a main contribution $H_0$ (e.g. the Coulomb potential) and a perturbation $H'$ which is small compared to $H_0$ and that will therefore just slightly change the main solution.
You then assume that you can expand the correct solution for the energy and the wavefunction as a perturbative series, i.e. in terms that are smaller and smaller: $$ |\psi_n\rangle = |\psi_n\rangle^0 + |\psi_n\rangle^1 + |\psi_n\rangle^2 ...$$
$$ E_n = E_n^0 + E_n^1 + E_n^2 ... $$
where the exponents signify the order of the term. Higher order terms are smaller, and therefore only needed in you want higher precision.
Now, the full TISE becomes: $$ (H_0 + H')\,(|\psi_n\rangle^0 + |\psi_n\rangle^1 + |\psi_n\rangle^2 ... ) = (E_n^0 + E_n^1 + E_n^2 + ...)\,(|\psi_n\rangle^0 + |\psi_n\rangle^1 + |\psi_n\rangle^2+...)$$
Now you take the $0^{th}$ order equation -- where each term is of $O^{th}$ order:
$$H_0 |\psi_n\rangle^0\rangle = E_n^0 |\psi_n\rangle^0,$$
which is just the unperturbed equation, and therefore the starting point for any perturbative calculation.
Now look at the $1^{st}$ order equation -- remember that two $1^{st}$ order terms multiplied together give you a $2^{nd}$ order terms, whereas you only want to keep the $1^{st}$ order ones. $H'$ is first order:
$$ H_0 |\psi_n\rangle^1 + H'|\psi_n\rangle^0 = E_n^0|\psi_n\rangle^1 + E_n^1 |\psi_n\rangle^0.$$
What you are after is $E_n^1$, i.e. the $1^{st}$ order contribution to the energy.
Multiply by $^0\langle \psi_n|$ from the left:
$$ ^0\langle \psi_n|H_0 |\psi_n\rangle^1 + ^0\langle \psi_n|H'|\psi_n\rangle^0 = ^0\langle \psi_n|E_n^0|\psi_n\rangle^1 + ^0\langle \psi_n|E_n^1 |\psi_n\rangle^0,$$
$$E_n^0 \,^0\langle \psi_n|\psi_n\rangle^1 + \, ^0\langle \psi_n|H'|\psi_n\rangle^0 = E_n^0\,^0\langle \psi_n|\psi_n\rangle^1 + E_n^1 \,^0\langle \psi_n|\psi_n\rangle^0,$$
$$ \implies (E_n^0 - E_n^0) \,^0\langle \psi_n|\psi_n\rangle^1 + ^0\langle \psi_n|H'|\psi_n\rangle^0 = E_n^1 \,^0\langle \psi_n|\psi_n\rangle^0.$$
Assuming normalised states, $^0\langle \psi_n|\psi_n\rangle^0 = 1$, so:
$$ E_n^1 = ^0\langle \psi_n|H'|\psi_n\rangle^0 $$.
As other have noted, there's a mistake in you formula.
The same procedure applies for higher order corrections.
Best Answer
Spin1/2 particle
Ususally, in this kind of Hamiltonian, people uses $s=s_z$, where
$$s=s_z=\left[ \begin{array}{cc} 1 & 0 \\ 0 & -1\end{array} \right].$$
Then, your unperturbed hamiltonian $H_0$ is: $$H_0=-\mu s\cdot B_0 = -\mu \left[ \begin{array}{cc} 1 & 0 \\ 0 & -1\end{array} \right]B_{0,z}. $$
Then the eigen vectors of energy are:
$$|\psi^0_+\rangle=\left[ \begin{array}{c} 1 \\ 0\end{array} \right],$$ $$|\psi^0_-\rangle=\left[ \begin{array}{c} 0\\ 1 \end{array} \right].$$
Perturbation solution
Then you want to compute $|\psi_+\rangle$ and $|\psi_-\rangle$ for the perturbed Hamiltonian $H=H_0-\mu B_1 s_x$, where $$s_x=\left[ \begin{array}{cc} 0 & 1 \\ 1 & 0\end{array} \right].$$
As you said, you have to compute the following quantities (note I use $+,-$ instead of $n=0,1$. Which became:
$$\psi^{(1)}_+=\sum_{n\neq +} \psi^{(0)}_{n'}\frac{\langle\psi_{n'}^{(0)}|-\mu B_1s_x|\psi_{+}^{(0)}\rangle}{E_+^{(0)}-E_{n'}^{(0)}}=\psi^{(0)}_{-}\frac{\langle\psi_{-}^{(0)}|-\mu B_1s_x|\psi_{+}^{(0)}\rangle}{E_+^{(0)}-E_{-}^{(0)}}$$
$$\psi^{(1)}_-=\sum_{n\neq -} \psi^{(0)}_{n'}\frac{\langle\psi_{n'}^{(0)}|-\mu B_1s_x|\psi_{-}^{(0)}\rangle}{E_-^{(0)}-E_{n'}^{(0)}}=\psi^{(0)}_{+}\frac{\langle\psi_{+}^{(0)}|-\mu B_1s_x|\psi_{-}^{(0)}\rangle}{E_-^{(0)}-E_{+}^{(0)}}$$
Put here vectors and matrices we just found and let me know if you get zero.