Here is a possibility to work.

Let $A$ be a ring such that $2$ is invertible, for example $\Bbb Z/3$ or $\Bbb Z/9$. We denote by $\frac 12$ its inverse.
Assume $2$ is not a square in $A$. (This is the case in the above examples.)

Consider $A[a]$ to be the ring $A[Y]/(Y^2-2)$, where $A[Y]$ is the polynomial ring in $Y$ over $A$. (And $a$ is $Y$ modulo the ideal generated by $(Y^2-2)$. So formally, $a=\sqrt2$.

Now observe that the equation satisfied by the $x$ element in the OP, $(x+a)^2=2(x+a)$, can be written as
$$(x+a)(x+a-2)=0\ .$$
The two factors are relatively prime, in fact
$$\frac 12(x+a)-\frac 12(x+a-2) = 1\ .$$
This means that the ring
$$
R=A[a][X]\ /\ (\ (X+a)(X+a-2)\ )
$$
is isomorphic to two copies of $A[a]$ via the map
$$
R\to A[a]\times A[a]\ ,\ f(X)\to(f(-a), f(2-a))\ .
$$
The inverse map takes $(s,t)\in A[a]\times A[a]$ and maps it into $\frac 12t(X+a)-\frac 12s(X+a-2)$.

Now we have to solve $Z^2=1$ in the ring $A[a]\times A[a]$. This can be done also manually. An element of the shape $Z=u+va\in A[a]$, $u,v\in A$ satisfies $Z^2 = 1$ iff $(u^2+2v^2)+2auv=1$. In case $A$ has zero divisors we have to take care of $uv=0$ somehow. The possible values for $u,v$ that may lead to a solution satisfy at any rate $u^3=u$ and $2v^3=v$. Together with $Z$, also its "conjugate" $\bar Z=u-va$ is also a solution, and the "norm" $N(Z)=Z\bar Z=(u+va)(u-va)=u^2-2v^2$ is $1$. So it is a good idea to search for solutions of this "Pell equation" over $A$.

Let us now write some lines of code for the given case.

The brute force search is:

```
r = Zmod(9)
R.<a,x> = PolynomialRing(r, order='lex')
J = R.ideal( [a^2-2, (x+a)^2-2*(x+a)] )
S = R.quotient(J)
for [s, t, u] in cartesian_product([r, r, r]):
Z = S(s + t*a + u*x)
if Z^2 == S(1):
print(f"Z = {s} + {t} a + {u} x")
```

Results:

```
Z = 1 + 0 a + 0 x
Z = 1 + 8 a + 8 x
Z = 8 + 0 a + 0 x
Z = 8 + 1 a + 1 x
```

This fits with the situation of finding all points $Z=(Z_1,Z_2)$ over the ring $R=A[a]=\Bbb Z/9[a]=\Bbb Z/9[\sqrt 2]$ with $Z^2=(1,1)=1_R$.

```
sage: r = Integers(9)
sage: R.<Y> = r[]
sage: Q.<a> = R.quotient(Y^2-2)
sage: a^2
2
sage: for r1, r2 in cartesian_product([r, r]):
....: Z1 = r1 + r2*a
....: if Z1^2 == Q(1):
....: print(Z1)
....:
1
8
```

These are the only Hensel lifts of the corresponding units over $\Bbb Z/3$:

```
sage: U.<A> = PolynomialRing(GF(3))
sage: F.<a> = GF(3^2, modulus = A^2-2)
sage: a^2
2
sage: [ f for f in F if f^2 == F(1) ]
[2, 1]
sage: # of course, this is a field...
```

Just make a triple loop.

@FrédéricC Yes, in this particular case I can do it, but I'd like to learn how to tackle this kind of computation more in general. I'm going to need computations over infinite rings too. The specific mention to $\mathbb{Z}_9$ is in case that some module does not work well over rings which are not fields.

@FrédéricC In fact that presentation of the problem was a mistake on my part, since I don't know if the linear combinations are going to span the whole ring. I need some function giving the set of solutions reduced by the Gröbner basis. Thank you for your help!