It might be useful to write out all entries of those matrices; while your notation generalises nicely to higher dimensions, it's at first confusing that the $b_i, c_i$ etc. stand for various matrices of various different shapes. Written entirely out, a general element of that split form of type $B_2$, often called $\mathfrak{so}(3,2)$, looks like this:
$$\pmatrix{0&e&f&g&h\\
-g&a&c&0&i\\
-h&d&b&-i&0\\
-e&0&j&-a&-d\\
-f&-j&0&-c&-b}$$
where all ten variables are from the ground field.
The diagonal matrices contained in this form a Cartan subalgebra, and let's say we call $\alpha$ the root sending
$$\alpha: \pmatrix{0&0&0&0&0\\
0&a&0&0&0\\
0&0&b&0&0\\
0&0&0&a&0\\
0&0&0&0&b} \mapsto a-b ;$$
and $\beta$ the one sending the above element to $-a$; the other roots are $\alpha+\beta$, $\alpha+2\beta$, and all their negatives. E.g. the root space to $\beta$ is the one where everything except the letter $e$ is zero; the root space to $\alpha+\beta$ is the one where everything except $h$ is zero; the root space to $\alpha+2\beta$ consists of the matrices where everything but $j$ is zero, etc.
Let's believe the relations in Dietrich Burde's answer. How to find an explicit set of corresponding matrices?
Well, the $h$'s in any Cartan-Weyl basis are already determined by $\rho(h_\rho)=2$ and living in $[\mathfrak g_\alpha, \mathfrak g_{-\alpha}]$, giving us here
$$h_\alpha= \pmatrix{0&0&0&0&0\\
0&1&0&0&0\\
0&0&-1&0&0\\
0&0&0&-1&0\\
0&0&0&0&1},
h_\beta= \pmatrix{0&0&0&0&0\\
0&-2&0&0&0\\
0&0&0&0&0\\
0&0&0&2&0\\
0&0&0&0&0}.$$
Further, $x_\beta$ is some
$$x_\beta= \pmatrix{0&x&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0\\
-x&0&0&0&0\\
0&0&0&0&0}$$
with $x \neq 0$ in the ground field; because we want $[x_{-\beta}, x_\beta]=h_\beta$, we need
$$x_{-\beta}= \pmatrix{0&0&0&2/x&0\\
-2/x&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0}.$$
Likewise but easier we have
$$x_\alpha= \pmatrix{0&0&0&0&0\\
0&0&y&0&0\\
0&0&0&0&0\\
0&0&0&0&0\\
0&0&0&-y&0}, x_{-\alpha}= \pmatrix{0&0&0&0&0\\
0&0&0&0&0\\
0&-1/y&0&0&0\\
0&0&0&0&1/y\\
0&0&0&0&0}$$
for some $y\neq 0$. Now $x_{\alpha+\beta}$ can be determined as
$$x_{\alpha+\beta} = [x_\beta, x_\alpha]= \pmatrix{0&0&xy&0&0\\
0&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0\\
-xy&0&0&0&0}$$
and $x_{\alpha+2\beta}$ as
$$x_{\alpha+2\beta} = \frac12 [x_\beta, x_{\alpha+\beta}]= \pmatrix{0&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0\\
0&0&-\frac{1}{2}x^2y&0&0\\
0&\frac{1}{2}x^2y&0&0&0}.$$
To find a good $x_{-\alpha-\beta}$, first set a general candidate for it as
$$\pmatrix{0&0&0&0&z\\
0&0&0&0&0\\
-z&0&0&0&0\\
0&0&0&0&0\\
0&0&0&0&0}.$$
But we need $[x_\beta, x_{-\alpha-\beta}]=-2x_{-\alpha}$ and/or $[x_\alpha, x_{-\alpha-\beta}]=x_{-\beta}$, both forcing $z=\frac{2}{xy}$.
I stop here because I'm confident we could continue and find the last few $x_\rho$'s; we end up with two arbitrary parameters $x,y$, which of course we could ocnveniently set to $1$.