First of all, we assume that we have a finite MDP, i.e. the set of states $\mathcal{S}$, the set of actions $\mathcal{A}$ and the set of rewards $\mathcal{R}$ all have a finite number of elements (I didn't think about how the explanations below would extend to other cases, but I suspect you will need differential equations).
For simplicity, let's only consider the value function $v$ (as opposed to the state-action value function $q(s, a)$, but this also applies to $q$). The value function $v$ is defined for all states, i.e. it's a function of the form $v : \mathcal{S} \rightarrow \mathbb{R}$ or, with an alternative notation, $v(s), \forall s \in \mathcal{S}$. So, we can define this function as a vector $\mathbf{v}$ of dimension $|\mathcal{S}| = n$, i.e. $\mathbf{v} \in \mathbb{R}^{|\mathcal{S}|}$, where the $i$th element contains the value of the $i$th state (so we need a function that maps states to indices of this vector, but this is trivial).
The fact that you can represent the value function, in this finite MDP, as a vector should already suggest that you can find this value function by solving a linear system of equations.
However, let me show that by starting with the definition of the value function you also provided
\begin{align}
v_{\pi}(s)
&=
\sum_{a} \pi(a\rvert s) \sum_{s^{\prime}}\sum_{r} p(s^{\prime}, r \rvert s,a)[r + \gamma v_{\pi}(s^{\prime})] \label{1}\tag{1}, \; \forall s \in \mathcal{S}
\end{align}
which can be expanded as follows
\begin{align}
v_{\pi}(s)
&=
\sum_{a} \pi(a\rvert s) \sum_{s^{\prime}}\sum_{r} \left[p(s^{\prime}, r \rvert s,a)r + \gamma p(s^{\prime}, r \rvert s,a) v_{\pi}(s^{\prime}) \right] \\
&=
\sum_{a} \pi(a\rvert s) \left[ \sum_{r} \underbrace{\sum_{s^{\prime}} p(s^{\prime}, r \rvert s,a)}_{\text{Marginalization of }p \text{ over } s'}r + \gamma \sum_{s^{\prime}} \underbrace{\sum_{r} p(s^{\prime}, r \rvert s,a) }_{\text{Marginalization of }p \text{ over } r} v_{\pi}(s^{\prime}) \right]\\
&=
\sum_{a} \pi(a\rvert s) \left[ \sum_{r} p(r \rvert s, a)r + \gamma \sum_{s^{\prime}}p(s^{\prime} \rvert s,a) v_{\pi}(s^{\prime}) \right] \\
&=
\sum_{a} \pi(a\rvert s) \left[ r(s, a) + \gamma \sum_{s^{\prime}}p(s^{\prime} \rvert s,a) v_{\pi}(s^{\prime}) \right] \label{2}\tag{2}, \; \forall s \in \mathcal{S}
\end{align}
where
- $\sum_{r} p(r \rvert s, a)r = r(s, a)$ (see this).
- $\sum_{s^{\prime}}p(s^{\prime}, r \rvert s,a) = p(r \rvert s, a)$ (marginalization)
- $\sum_{r} p(s^{\prime}, r \rvert s,a) =p(s^{\prime} \rvert s,a) $ (marginalization)
In this form, as in equation \ref{2}, the value function can also be written in a different notation
\begin{align}
v_{\pi}(s)
&=
\sum_{a} \pi(a\rvert s) \left[ R_{s}^a + \gamma \sum_{s^{\prime}}P_{ss'}^a v_{\pi}(s^{\prime}) \right] \label{3}\tag{3}, \; \forall s \in \mathcal{S}
\end{align}
where
- $R_{s}^a = r(s, a)$
- $P_{ss'}^a = p(s^{\prime} \rvert s,a)$
We can still write equation \ref{3} in a "simpler" form as follows
\begin{align}
v_{\pi}(s)
&=
\sum_{a} \pi(a\rvert s) R_{s}^a + \gamma \sum_{s^{\prime}} \sum_{a} \pi(a\rvert s) P_{ss'}^a v_{\pi}(s^{\prime}) \\
&=
R_{s}^\pi + \gamma \sum_{s^{\prime}} P_{ss'}^\pi v_{\pi}(s^{\prime})
\label{4}\tag{4}, \; \forall s \in \mathcal{S}
\end{align}
where
- $\sum_{a} \pi(a\rvert s) R_{s}^a = R_{s}^\pi$
- $\sum_{a} \pi(a\rvert s) P_{ss'}^a = P_{ss'}^\pi$
We can write the definition of the value function in \ref{4} in matrix form for all states $s \in \mathcal{S}$ as follows
\begin{align}
\begin{bmatrix}
v_\pi(1) \\
\vdots \\
v_\pi(n)
\end{bmatrix}=
\begin{bmatrix}
{R}_1^\pi \\
\vdots \\
{R}_n^\pi
\end{bmatrix}
+\gamma
\begin{bmatrix}
{P}_{11}^\pi & \dots & {P}_{1n}^\pi\\
\vdots & \ddots & \vdots\\
{P}_{n1}^\pi & \dots & {P}_{nn}^\pi
\end{bmatrix}
\begin{bmatrix}
v_\pi(1) \\
\vdots \\
v_\pi(n)
\end{bmatrix}
\tag{5}\label{5},
\end{align}
which can be written in a more compact form as follows
\begin{align}
\mathbf{v} = \mathbf{r} + \gamma \mathbf{P}\mathbf{v} \tag{6}\label{6},
\end{align}
which is a very compact form of the Bellman equation (which is a recursive equation: as you can notice, the $\mathbf{v}$ appears on the left and right of the equals sign) that represents the value function (i.e. the value function can be defined as a recursive equation).
In equation \ref{5}, the unknowns are the $|\mathcal{S}| = n$ values of the value function $v$ and there are $n$ equations, so it should now be clear why we can solve this problem by solving a system of equations. Note that here it's assumed that $\pi$, $r(s, a)$ and $p$ are given and known, which, generally, is not the case, that's why we use algorithms like Q-learning.
