Solving equations without side effects
As far as I know, I think it is possible to build a functional implementation of solve that does behave as a deterministic function, with no side effects.  Yes, there might be multiple possible outputs that are all valid.  However, I think it's possible to pick an implementation where it always outputs the same answer.  (Certainly, the existence of multiple possible outputs doesn't on its own render it impossible to build such an implementation.)

Here is a simpler example, to give you some intuition.  Suppose we want to find a solution to $|x|=c$, where $c$ is a constant given in advance.  Then there are two solutions ($x=+c$ and $x=-c$).  Nonetheless, you can still build a functional program to output a solution to this, i.e., on input $c$, the code returns a solution to $|x|=c$.  An example implementation is the identity function (always return whatever the input was).  This is functional and has no side effects.

So just because there are multiple valid outputs doesn't require the code to randomly or nondeterministically choose among them.  Often you can write code that outputs one of them, and is consistent and deterministic about which one it outputs.

For instance, here is a plausible such implementation for your particular problem:


Solve the linear system "minimize $x_1$ such that $Ax=b$" using linear programming.  (Here I write $x_1$ for the first variable of $x$, i.e., I am assuming $x=(x_1,\dots,x_n)$.)  (If it has no solution, return "No solution".)  Let $u^{(1)}$ denote the resulting such solution.
Solve the linear system "minimize $x_2$ such that $Ax=b$ and $x_1 = u^{(1)}_1$" using linear programming.  Let $u^{(2)}$ denote the resulting such solution.
Solve the linear system "minimize $x_3$ such that $Ax=b$ and $x_1 = u^{(1)}_1$ and $x_2 = u^{(2)}_2$" using linear programming.  Let $u^{(3)}$ denote the resulting such solution.
Repeat, until you have obtained $u^{(n)}$, where $n$ is the number of variables.  Return $u^{(n)}$.


Note that this algorithm returns the lexicographically least solution to $Ax=b$.  Because linear programming is in P, I believe each "solve the linear system" step can be implemented with a deterministic, functional algorithm.