I was able to reverse-engineer the method of Matsunaga Ryohitsu and find that it does work, modulo understanding the following two points about it, that are not necessarily clear from the wording (either Mikami's or Dickson's):
When it is impossible to subtract a $b_i$ from the current remainder $A$, what we must add to $A$ is the next unused $a_j$ (after the latest $a_{j-1}$ used), which is not necessarily the "corresponding" $a_i$ (i.e., $j$ may be a smaller index than $i$). Alternatively, instead of changing $a_{j+1} = a_j - 2$ at every step, we change it only when required. That is, the $a_i$'s and $b_i$'s aren't meant to be recalculated simultaneously at each step, but the next $a_i$ is only to be picked when necessary. This is clearer from Mikami's wording than Dickson's.
Mikami's wording "let the values of $a$ and $b$ employed in the last place be $a'$ and $b'$" (or Dickson's "$a', b'$ are the values last employed") are probably not what the original intended, which seems to be to take $a'$ and $b'$ to be the next ones after the ones last used in an addition or subtraction. This can be sort of coerced to fit the original wording, if we imagine that as soon as we use up an $a_j$ or a $b_i$, we write the next ones immediately, and so $a'$ and $b'$ are taken to be the last ones writen down.
With that in mind, the method can be viewed as follows. We want $x^2 + y^2 = k$, so we want to search over all pairs of positive integers $(x, y)$ until we find a pair satisfying the equation. We will organize the search as follows. Let $\frac{k}2 = r^2 + R$ (where $r^2$ is the largest square not greater than $\frac12k$, and $R$ is the "remainder") so that $k = 2r^2 + 2R$. If we assume wlog that $x \le y$, we see that we must have $\color{red}{y \ge r}$ (because if $x \le y < r$, then $x^2 + y^2 < 2r^2 \le k$) and that $\color{red}{x \le r}$ (because if $r < x \le y$, then $x^2 + y^2 \ge 2(r+1)^2 > k$).
So let's start with $x, y$ both equal to $r$, for which $x^2 + y^2 = 2r^2 \le k$.
We will successively increment $y$ by $1$ at a time as long as $x^2 + y^2 \le k$, then when we have overshot $k$ we will decrement $x$ by $1$ and continue, and so on. It is clear that this method will always work and find a solution if one exists, as it's an exhaustive search over all potential solutions $(x,y)$: for every $y \ge r$, we try the largest $x$ such that $x^2 + y^2 \le k$.
The beauty of this method is that accomplishes the exhaustive search without any squaring or multiplications, just additions and subtractions, that too of relatively smaller numbers.
More concretely, for a fixed $x$, suppose we increase $y$ by $1$. When we replace $y$ with $y+1$, the left-hand side (namely $x^2 + y^2$) increases to $x^2 + (y+1)^2$, i.e., it increases by $2y + 1$: twice the current candidate value for $y$, plus $1$. Or, if we consider candidates for $2y$ instead, the left-hand-side increases by the current candidate (for $2y$) plus $2$. This is what the $b$'s are doing in the method. (As the final $y = \frac12(b' - 1)$, this means $b' = 2y + 1$.)
Similarly, when we increase $y$ by too much, the left-hand side $x^2 + y^2$ will overshoot the right-hand side, and to compensate, it is time for try the next lower value for $x$ (or equivalently, for $2x$). When we replace $x$ with $x-1$, the left-hand side changes from $x^2 + y^2$ to $(x-1)^2 + y^2$, decreasing by $2x-1$, or twice the current value minus $1$. Or, if we consider candidate values for $2x$, it decreases by the current value for $2x$ minus $2$. These are the $a_j$'s, shifted by $1$. (As the final $x = \frac12(a' + 1)$, this means that $a' = 2x - 1$.)
A complete description of the method
We maintain, at all times $t$, a pair of integers $(x_t, y_t)$ such that $x_t^2 + y_t^2 \le k$. However, we won't work directly with the numbers $(x_t, y_t, k)$ except to output a solution; instead we work with the numbers
$$\begin{align}
a_t &= 2x_t - 1 ,\\
b_t &= 2y_t + 1, \\
A_t &= k - (x_t^2 + y_t^2).
\end{align}$$
Initially, at $t = 1$, (with the idea that $x_t = y_t = r$) set $a_t = 2r - 1$, $b_t = 2r + 1$, and $A_t = k - 2r^2$, where $r^2$ is the largest square not greater than $\frac{k}2$.
If $A_t = 0$, output the solution $x_t = \frac12(a_t + 1)$, $y_t = \frac12(b_t - 1)$.
While $A_t < b_t$, [need to decrease $x_t$ to make room for $y_t + 1$]
- Set $A_{t+1} = A_t + a_t$ and $a_{t+1} = a_t - 2$, keeping $b_{t+1} = b_t$.
- Increase $t$ to $t+1$.
[Now that $A_t \ge b_t$, can update $y_{t+1} = y_t + 1$.]
Set $A_{t+1} = A_t - b_t$ and $b_{t+1} = b_t + 2$, keeping $a_{t+1} = a_t$.
Update $t$ to $t+1$ and go back to step 2.
(Note that when we output a solution, we don't have to stop. If we keep going, we'll find all other solutions as well, if they exist. We can stop when $a_t < 0$.)
Some examples of the method (we don't have to restrict ourselves to even $k$), including all the ones in your question and more:
$k = 5$. In this case, $\frac12k = \frac52 = 1^2 + \frac32$, so $r = 1$ and $A_1 = 2R = 3$. Start with $a_1 = 2r - 1 = 1$, and $b_1 = 2r + 1 = 3$. Subtract, from $A = 3$, the current $b$: you already get $0$, so $a' = a_1 = 1$, and $b' = b_2 = b_1 + 2 = 5$. Now $x = \frac12 (a' + 1) = 1$ and $y = \frac12 (b' - 1) = 2$ is indeed a solution.
$k = 10$. In this case, $\frac12k = 5 = 2^2 + 1$, so $r = 2$ and $2R = 2$. Start with $a_1 = 2r - 1 = 3$ and $b_1 = 2r + 1 = 5$. It is already impossible to subtract $b_1$ from $A = 2R$, so add $a_1$ and carry out the subtraction: we get the new remainder $A$ to be $2R + a_1 - b_1 = 2 + 3 - 5 = 0$. So $a' = a_2 = a_1 -2 = 1$, and $b' = b_2 = b_1 + 2 = 7$. Now $x = \frac12(a' + 1) = 1$ and $y = \frac12(b' - 1) = 3$ is indeed a solution.
$k = 13$. We have $\frac12k = 2^2 + \frac52$, so $r = 2$ and $2R = 5$. $a_1 = 2r - 1 = 3$, $b_1 = 5$. Subtracing $b_1$ from $A=2R$ gives remainder $0$, so $a' = a_1 = 3$, and $b' = b_2 = 7$. This gives the solution $x' = \frac12(a' + 1) = 2$ and $y' = \frac12(b' - 1) = 3$.
$k = 17$. We have $\frac12k = 2^2 + \frac92$, so $r = 2$ and $2R = 9$. $a_1 = 2r - 1 = 3$, $b_1 = 5$. Update $A=2R$ to $A - b_1 = 4$. Next we cannot subtract $b_2 = 7$ from it, so add $a_1 = 3$ and subtract. Remainder is $0$, so $a' = a_2 = 1$, $b' = b_3 = 9$. This gives the solution $x' = \frac12(a' + 1) = 1$ and $y' = \frac12(b' - 1) = 4$.
$k = 41$. In this case, $\frac12k = 4^2 + \frac92$, so $r = 4$ and $2R = 9$. Start with $a_1 = 2r - 1 = 7$, and $b_1 = 2r + 1 = 9$. Subtracting $b_1$ from $A = 2R$ already gives remainder $0$, so $a' = a_1 = 7$ and $b' = b_2 = b_1 + 2 = 11$. Now $x' = \frac12(a' + 1) = 4$ and $y' = \frac12(b' - 1) = 5$ is indeed a solution.
For something less trivial, consider $k = 218$. In this case, $\frac12k = 109 = 10^2 + 9$, so $r = 10$ and $2R = 18$. Start with $a_1 = 2r - 1 = 19$, and $b_1 = 2r + 1 = 21$. Subtracting $b_1$ from $A = 2R$ is already impossible, so add $a_1$ and carry out the subtraction: the new value of $A$ is $2R + a_1 - b_1 = 16$, meanwhile we compute new values $a_2 = a_1 - 2 = 17$ and $b_2 = b_1 + 2 = 23$. Again, subtracting $b_2$ from $A$ is impossible, so we add $a_2$, thus setting $A$ now to $16 + a_2 - b_2 = 10$. Next compute $a_3 = a_2 - 2 = 15$, and $b_3 = b_2 + 2 = 25$. This time, subtracting $b_3$ is still not possible, but magically, it works if we add $a_3$, and indeed we get remainder $0$. So $a' = a_4 = a_3 -2 = 13$, and $b' = b_4 = b_3 + 2 = 27$. Now, $x' = \frac12(a' + 1) = 7$ and $y' = \frac12(b' - 1) = 13$ indeed a solution.
