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.