For an $n$ qubit system with $k$ linearly independent stabilizers, you can write out the $k\times 2n$ binary matrix that specifies the stabilizers: each row is one of the stabilizers. The first $n$ entries specify where $X$s are applied. The second $n$ entries specify where the $Z$s are applied. So, in your example case, you have
$$
H=\begin{bmatrix}
1 & 0 & 0 & 1 & 0 & 0 & 1 & 1 & 0 & 0 \\
0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1 & 0 \\
0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 1 & 0 \\
0 & 0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1
\end{bmatrix}.
$$
You're trying to find en error $e$ that satisfies
$$
He=s.
$$
So, what you want to do is use row operations (modulo 2) to cast $H$ into a standard form such that you can identify an entire $m\times m$ identity block (they don't have to be consecutive columns). You can always do this because the rows are linearly independent. You'll want to keep track of these row operations in a matrix $M$, i.e. $MH$ is cast into the standard form. Thus, now you're trying to solve
$$
MHe=Ms.
$$
The trick now is that the columns of $H$ corresponding to that identity block, you can just set the corresponding bits of $e$ equal to $Ms$.
In fact, in your example, you don't need an $M$; it's already in standard form. You're just looking at columns 1, 6, 5, 10. So, for any $s$, you're just going to set
$$
\begin{bmatrix}
e_1 \\
e_6 \\
e_5 \\
e_{10}
\end{bmatrix}=s.
$$
Thus, your error operator will look like
$$
Z_1^{s_1}Z_5^{s_3}X_1^{s_2}X_5^{s_4}
$$
(because in the error term, $e$, you've swapped the $X$ and $Z$ terms because $X$ stabilizers detect $Z$ errors).
In terms of coding, I often use Mathematica for this stuff because many of its functions have a built in Modulus 2 option. In python, however, you might consider the ldpc package. This has a mod2 sub-package that can do the inversion etc for you. Or, when it performs row reduction for you, it will also output the matrix $M$ that does the row reduction, which many packages don't.
