To convert this into a generic linear algebra problem you'd rewrite it in the form $Ax=b$ where $A$ is a matrix of stoichiometric coefficients of size $m \times n$; $x$ is a vector of length n of unknown fitting parameters (the heats you want to determine) and $b$ a known vector of length n.
For the example problem the set of equations can be written in matrix form as
$$\left[\begin{array}{rrrrr} 0 & -4 & -3 & 6 & 2\\ -2 & 0 & -1 & 2 & 0 \end{array}\right]\left[ \begin{array}{l} \Delta_\mathrm{f} H(\ce{H_2(g)})\\\Delta_\mathrm{f} H(\ce{NH_3(g)})\\\Delta_\mathrm{f} H(\ce{O_2(g)})\\\Delta_\mathrm{f} H(\ce{H_2O(l)})\\\Delta_\mathrm{f} H(\ce{N_2(g)}) \\ \end{array}\right] = \left[ \begin{array}{c } \Delta _\mathrm{r}H_1\\\Delta_\mathrm{r}H_2 \end{array}\right]$$
Since many of the terms in the vector of heats on the left-hand-side are known and equal to zero this can be rewritten as
$$\left[\begin{array}{rrlll} -4 & 6 \\ 0 & 2 \end{array}\right]\left[ \begin{array}{c } \Delta_\mathrm{f} H(\ce{NH_3(g)})\\\Delta_\mathrm{f} H(\ce{H_2O(l)}) \end{array}\right] = \left[ \begin{array}{c } \Delta_\mathrm{r} H_1\\\Delta_\mathrm{r} H_2 \end{array}\right]$$
(alternatively the original equation can be written as a $5 \times 5$ array by including additional trivial equations defining the heats of formation of the diatomic gases).
Solution of the above system should be easy, for instance in MATLAB:
A = [-4 6; 0 2];
b = [-1516 -572]';
x = A\b
The problem can of course be solved by inspection. The second equation can be solved for $\Delta_\mathrm{f} H(\ce{H_2O(l)})$, which can be inserted into the first to solve for $\Delta_\mathrm{f} H(\ce{NH3(g)})$.