Here is a brief sketch of how one can solve linear systems modulo any fixed integer$~n>0$. If $n$ is prime, usual Gaussian reduction of the system works. If $n$ is a product of distinct relatively prime factors$~n_i$, then by the Chinese remainder theorem the system is equivalent to the collection of systems obtained from it by reducing modulo$~n_i$ for every $i$. Supposing each of these systems can be solved (or proven inconsistent, which I consider a form of solving) one can combine the solutions to get the solution for the original system via the Chinese remainder theorem. In detail: if the system reduced modulo any$~n_i$ is inconsistent, so it the original system; otherwise there are parametrised solutions modulo each$~n_i$, one may need to introduce dummy parameters (occurring with coefficient$~0$ everywhere) in some of them to force the same number of free parameters for each$~n_i$, and then the CRT can be applied to a particular solution and to each of the parameters to lift each of them to a value modulo$~n$, which gives a parametrised solution over $\def\Z{\Bbb Z}\Z/n\Z$.
One can thus reduce to the case where each $n_i$ is a prime power; one still needs to bridge the gap between prime powers $p^k$ and the prime$~p$. For this I suggest using a form of Hensel lifting* to the system. First reduce the system modulo$~p$ and solve that. Then if there are any solutions, lift one of them arbitrarily to a non-solution modulo$~p^2$. Then add an unknown multiple of$~p$ as correction term to it, and express the equation that the sum become a true solution modulo$~p^2$; since any defect of the lifted non-solution is a multiple of$~p$ one can divide a factor$~p$ from the equations obtained, and get a linear system modulo$~p$ that can again be solved using Gaussian elimination. Once lifted to$~p^2$, continue similarly to lift to higher powers of$~p$, until reaching$~p^k$.
There is a complication to this idea for which I do not have an answer right now. Solutions will in general not be unique but parametrised; in principle this means one must apply the lifting separately to a particular solution and to the parametrised part (the solution to the corresponding homogeneous system). However it looks like sometimes the question of whether the particular solution can be lifted at all to a higher power of$~p$ may depend on the choice of that particular solution. If this occurs, it would would mean that during the lifting process one would need to allow for adapting the particular solution to avoid inconsistencies in the set of equations obtained. This would be an additional headache, but maybe it can be shown that this never occurs. One thing that certainly can occur is that a solution exists modulo$~p$ but none exist modulo$~p^2$; this unlike the situation of Hensel's lemma (which only deals with one equation, but which is non-linear). A trivial example of that situation is the system of equations $(x\equiv1,x\equiv3)\pmod4$ which is clearly inconsistent but whose reduction modulo$~2$ is not.
*This is actually due to Schönemann, a relatively unknown 19th century mathematician who also first formulated "Eisenstein's criterion" for irreducibility of polynomials.
In modulo arithmetic, "division" means multiplying by the multiplicative inverse, e.g., $b = \frac{1}{a}$ means the value which when multiplied by $a$ gives $1$ modulo the value, e.g., $ba \equiv 1 \pmod n$. Note you may sometimes see $b = a^{-1}$ instead to avoid using explicit "division". This works, and gives a unique value, in any cases where the value you're dividing and the modulus are relatively prime.
More generally, it'll work in all cases of $\frac{c}{a} \equiv e \pmod n$ where $d = \gcd(a,n)$ and $d \mid c$ since this gcd value "cancels" in the division. Thus, the resulting equivalent modulo equation of $\frac{c'}{a'} \equiv e \pmod n$, where $c' = \frac{c}{d}$ and $a' = \frac{a}{d}$ has $\gcd(a',n) = 1$, has a solution. However, as Bill Dubuque's comment indicates, this assumes you're doing integer division to the extent of removing the common factor of $d$. Note that $a^{-1}$ doesn't exist modulo $n$ in this case. However, $(a')^{-1}$ does exist modulo $\frac{n}{d}$, so a possible interpretation would be $\frac{c'}{a'} \equiv c'(a')^{-1} \equiv e \pmod{\frac{n}{d}}$.
As for why the multiplicative inverse $b = a^{-1}$ exists modulo $n$ if $\gcd(a,n) = 1$, Bézout's identity states that in such cases there exist integers $x$ and $y$ such that
$$ax + ny = 1 \tag{1}\label{eq1}$$
As such $ax \equiv 1 \pmod n$ so $x \equiv a^{-1} = b \pmod n$. This value must be unique, modulo $n$, because if there was another value $x'$ such that $xa \equiv x'a \equiv 1 \pmod n$, then $(x - x')a \equiv 0 \pmod n$. Since $\gcd(a,n) = 1$, this means that $x - x' \equiv 0 \pmod n \; \iff \; x \equiv x' \pmod n$.
Bézout's identity also shows that if $a$ and $n$ are not relatively prime, e.g., $\gcd(a,n) = d \gt 1$, then \eqref{eq1} becomes
$$ax + ny = d \tag{2}\label{eq2}$$
with the integers of the form $ax + ny$ are always multiples of $d$, so it can't be congruent to $1$ and, thus, $a$ would not have a multiplicative inverse.
Best Answer
There's also a "cheats" method available here. There are only $25$ possible values of $(x,y) \in (\mathbb{Z}_5)^2$. We can just check them one-by-one, and see which ones work.
We could do this by hand, or on a computer. In GAP:
returns the single solution $(x,y)=(2,4)$.