I have worked through two rather explicit approaches to computing the
Massey product. One way is to follow the notes of Bruner, using
theorem 2.6.5 and the paragraphs which follow. The notes are very good, and there isn't much I could possibly add to them. The second is to use
the cobar complex and work out the Massey product explicitly. For this
approach, the computation is much simpler if we use the results in
McCleary's book "A User's Guide to Spectral Sequences," see pages
378--379, and 128--130.
Let $\mathcal{A}$ denote the Steenrod algebra at the prime 2. Recall
$\mathcal{A}$ is generated by the Steenrod squares $Sq^i$ for $i\geq
0$ and the admissible monomials form an $\mathbb{F}_2$-basis for
$\mathcal{A}$. I will use $-^{\vee}$ for the functor of graded
$\mathcal{A}$-modules $Hom^{*}_{\mathcal{A}}(-, \mathbb{F}_2)$. Recall
that $\mathcal{A}$ has an augmentation $\epsilon : \mathcal{A} \to
\mathbb{F}_2$. Denote $\ker(\epsilon)$ by $I(\mathcal{A})$.
We first need the (reduced) bar resolution of $\mathbb{F}_2$. This
takes the form
\begin{equation*}
\mathbb{F}_2 \leftarrow \mathcal{A} \leftarrow \mathcal{A} \otimes I(\mathcal{A}) \leftarrow \mathcal{A} \otimes I(\mathcal{A}) \otimes I(\mathcal{A}) \leftarrow \cdots
\end{equation*}
Let $B_n = \mathcal{A} \otimes I(\mathcal{A})^{\otimes n}$ to simplify
notation. The tensor products are over $\mathbb{F}_2$, and the
traditional notation for an element in $B_n$ is $\gamma[\gamma_1 \vert
\gamma_2 \vert \cdots \vert \gamma_n]\equiv \gamma \otimes \gamma_1
\otimes \cdots \otimes \gamma_n$. A reason for this is that we
consider $B_n$ to be an $\mathcal{A}$-module with action $x \cdot
\gamma[\gamma_1 \vert \cdots \vert \gamma_n] = (x\gamma)[\gamma_1
\vert \cdots \vert \gamma_n]$. The maps $d_n : B_n \to B_{n-1}$ in the
resolution are defined by
\begin{equation*}
d_n(\gamma[\gamma_1 \vert \cdots \vert \gamma_n]) = \gamma\gamma_1[\gamma_2 \vert \cdots \vert \gamma_n] + \gamma[\gamma_1\gamma_2 \vert \cdots \vert \gamma_n] + \cdots + \gamma[\gamma_1 \vert \cdots \vert \gamma_{n-2} \vert \gamma_{n-1}\gamma_{n}].
\end{equation*}
This is a resolution of $\mathbb{F}_2$ by free $\mathcal{A}$-modules,
as one can check. See pages 242--248 in McCleary for more
details. With this in hand, if we apply $-^{\vee} =
Hom^*_{\mathcal{A}}(-,\mathbb{F}_2)$ and take cohomology of the
resulting chain complex, we will get
$Ext^{**}_{\mathcal{A}}(\mathbb{F}_2,\mathbb{F}_2)$.
Let us now dualize the bar construction above, i.e., apply our functor
$-^{\vee}$. We will get
\begin{equation*}
\mathcal{A}^{\vee} \rightarrow (\mathcal{A} \otimes
I(\mathcal{A}))^{\vee} \rightarrow (\mathcal{A} \otimes I(\mathcal{A}) \otimes
I(\mathcal{A}))^{\vee} \rightarrow \cdots
\end{equation*}
which is
\begin{equation*}
\mathbb{F}_2 \rightarrow Hom^*_{\mathbb{F}_2}(I(\mathcal{A}),\mathbb{F}_2) \rightarrow Hom^*_{\mathbb{F}_2}(I(\mathcal{A}),\mathbb{F}_2) \otimes Hom^*_{\mathbb{F}_2}(I(\mathcal{A}),\mathbb{F}_2) \rightarrow \cdots.
\end{equation*}
An $\mathcal{A}$-module map $B_n \to \mathbb{F_2}$ is determined by its
values on an $\mathcal{A}$-basis, and the elements of the form
$1[\gamma_1 \vert \cdots \gamma_n]$ is one such basis, where all of
the $\gamma_i$'s are admissible monomials of positive degree. Thus we
see $B_n^{\vee} \cong Hom^*_{\mathbb{F}_2}(I(\mathcal{A}),
\mathbb{F}_2) = I(\mathcal{A})^{dual}$. Here we define $X^{dual} = Hom_{\mathbb{F}_2}(X, \mathbb{F}_2)$.
The dualized maps $d_n^{\vee} : B_{n-1}^{\vee} \to B_{}^{\vee}$ have a
very handy description if we use the description of the dual Steenrod
algebra $\mathcal{A}^{dual} = Hom_{\mathbb{F}_2}^*(\mathcal{A},
\mathbb{F}_2)$ determined by Milnor in his 1958 paper ``The Steenrod
algebra and its dual.'' Namely, $\mathcal{A}^{dual} =
\mathbb{F}_2[\zeta_1, \zeta_2, ... ]$ where $\mathrm{deg} \zeta_i = 2^i -
1$. Furthermore, the coproduct $\phi^{dual}$ for $\mathcal{A}^{dual}$
is determined by $\phi^{dual}(\zeta_i) = \sum_{j+k = i}
\zeta_j^{2^k}\otimes \zeta_k$. As $\phi^{dual}$ is an algebra
homomorphism, this determines the coproduct. To start off, $d_1^{\vee}
: \mathbb{F}_2 \to I(\mathcal{A})^{dual}$ is the zero map $1
\mapsto 0$. Furthermore, the map $d_2^{\vee} :
I(\mathcal{A})^{dual} \to I(\mathcal{A})^{dual}\otimes
I(\mathcal{A})^{dual}$ is given by $[x] \mapsto \sum [x' \vert x'']$ where
$\phi^{dual}(x) = 1\otimes x + x \otimes 1 + \sum x' \otimes x''$.
Here are some concrete computations. As $\phi^{dual}(\zeta_1) =
1\otimes \zeta_1 + \zeta_1 \otimes 1$, we have $\phi^{dual}(\zeta_1^2)
= \zeta_1^2 \otimes 1 + 1 \otimes \zeta_1^2$ and
\begin{equation*}
\phi^{dual}(\zeta_1^3) = 1\otimes \zeta_1^3 + \zeta_1 \otimes
\zeta_1^2 + \zeta_1^2 \otimes \zeta_1 + \zeta_1^3 \otimes 1
\end{equation*}
Therefore $d_2^{\vee}(\zeta_1) = 0$, $d_2^{\vee}(\zeta_1^2) = 0$ and
\begin{equation*}
d_2^{\vee}([\zeta_1^3]) = [\zeta_1 \vert \zeta_1^2] + [\zeta_1^2 \vert \zeta_1].
\end{equation*}
Another computation we will need is for $\zeta_2$. Recall $\zeta_2$ is
of degree $3$, and is the dual of $Sq^2Sq^1$. The coproduct of
$\zeta_2$ is $\phi^{dual}(\zeta_2) = 1\otimes \zeta_2 + \zeta_1^2
\otimes \zeta_1 + \zeta_2 \otimes 1$. Therefore
\begin{equation*}
d_2^{\vee}([\zeta_2]) = [\zeta_1^2 \vert \zeta_1].
\end{equation*}
Now what makes it possible to compute (and define I suppose!) the
Massey products in $Ext_{\mathcal{A}}(\mathbb{F}_2, \mathbb{F}_2)$ is
that the dual of the bar resolution---which is called the cobar
construction---is a DG algebra whose cohomology is
$Ext_{\mathcal{A}}(\mathbb{F}_2, \mathbb{F}_2)$. One product in the
bar resolution is the juxtaposition product (or concatenation
product). That is, define
\begin{equation*} [\gamma_1 \vert \cdots \vert \gamma_n] * [\beta_1
\vert \cdots \vert \beta_m] = [\gamma_1 \vert \cdots \vert \gamma_n
\vert \beta_1 \vert \cdots \vert \beta_m].
\end{equation*}
This turns out to agree with the composition product, which is a good
thing, see McCleary pg. 383. So with respect to the juxtaposition
product, the differentials $d_*^{\vee}$ satisfy the Leibniz rule. So
we can define the Massey products! The juxtaposition product is nice
because it allows us to do computations easily.
Let's now finally compute $\langle h_0, h_1, h_0 \rangle = \left\{[ a
*h_0 + h_0 *b] \, \vert \, d_2^{\vee}a = h_0*h_1 \, d_2^{\vee} b =
h_1 * h_0 \right\}$. In this definition, the brackets mean
cohomology class. We must first figure out what we mean by $h_0$ and
$h_1$. Let us compute $Ext_{\mathcal{A}}^{1,1}$ and
$Ext_{\mathcal{A}}^{1,2}$. For the first, the image is 0, and the
kernel of $d_2^{\vee}$ in degree 1 is $[\zeta_1]$. So let us define
$h_0 = [\zeta_1]$. Next, we see that the kernel of $d_2^{\vee}$ in
degree 2 is $[\zeta_1^2]$ and again the image is 0. So let us define
$h_1 = [\zeta_1^2]$. The product $h_0 * h_1$ is vanishes in
$Ext^{1,3}$, so it must be the boundary of something. By our above
calculations, we see that $d_2^{\vee}([\zeta_1^3] + [\zeta_2]) =
[\zeta_1 \vert \zeta_1^2]$. Also, the product $h_1 * h_0 = [\zeta_1^2
\vert \zeta_1]$ vanishes in $Ext^{1,3}$, so it must also be the
boundary of something. Our above calculations show that
$d_2^{\vee}([\zeta_2]) = [\zeta_1^2 \vert \zeta_1]$. We now can
identify one element in the Massey product. Namely:
\begin{equation*}
[\zeta_1^3 + \zeta_2 \vert \zeta_1] + [ \zeta_1 \vert \zeta_2] =
[\zeta_1^3 \vert \zeta_1] + [\zeta_2 \vert \zeta_1 ] + [\zeta_1 \vert \zeta_2].
\end{equation*}
Which element is this in $Ext^{2,4}$? Well, let us calculate some
more. First off,
\begin{align*}
\phi^{dual}(\zeta_2\zeta_1)& = (1\otimes \zeta_2 + \zeta_1^2 \otimes \zeta_1 + \zeta_2 \otimes 1)(1\otimes \zeta_1 + \zeta_1\otimes1)\\
& = 1\otimes \zeta_2\zeta_1 + \zeta_2\zeta_1 \otimes 1 + \zeta_1^3 \otimes \zeta_1
+ \zeta_2 \otimes \zeta_1 + \zeta_1^2 \otimes \zeta_2^2 + \zeta_2 \otimes \zeta_2,
\end{align*}
hence
\begin{equation*}
d_2^{\vee}([\zeta_2\zeta_1]) = [\zeta_1^3 \vert \zeta_1 ] + [\zeta_2 \vert \zeta_1] + [\zeta_1^2 \vert \zeta_1^2] + [\zeta_1\vert \zeta_2].
\end{equation*}
Thus we see that our class $[\zeta_1^3 \vert \zeta_1] + [\zeta_2 \vert
\zeta_1 ] + [\zeta_1 \vert \zeta_2]$ in $\langle h_0 , h_1,
h_0\rangle$ is cohomologous to $[\zeta_1^2 \vert \zeta_1^2 ] = h_1 *
h_1 = h_1^2$ which is a non-zero class in $Ext^{1,3}$. By further direct
computation, one can determine that the indeterminacy, which is
$h_0Ext^{2,3} + Ext^{2,3}h_0$, of the Massey product $\langle h_0, h_1,
h_0\rangle$ is $0$, so that we have completely determined it.