In your case, here are some results for you to either derive or prove through induction:
$$\sum_{i=1}^n 1 = n$$
$$\sum_{i=1}^n i = \frac{1}{2}n(n+1)$$
$$\sum_{i=1}^n i^2 = \frac{1}{6}n(n+1)(2 n+1)$$
$$\sum_{i=1}^n i^3 = \frac{1}{4}n^2(n+1)^2$$
Note that the summation is linear, so that for a term as you have above, you can break it up into
$$\sum_{i=1}^n 6 + \sum_{i=1}^n 2 i - \sum_{i=1}^n 2 = 6 \sum_{i=1}^n 1 + 2 \sum_{i=1}^n i - 2 \sum_{i=1}^n 1$$
I'll describe a rudimentary method for making the plots in Figure 4 of the paper. Since the three equations are so similar I'll just consider the first,
$$
-37.5 C/1000 - 373.41 B^{-0.6269} + 9.1256 \log(C+B) - 50.029 = 0. \tag{1}
$$
We'll solve $(1)$ for $C$ numerically as $B$ ranges over the interval $[600,3500]$. For each new application of Newton's method we'll use linear extrapolation to determine the new initial guess.
We need two starting points. Let's take $B = 600$ and define the function
$$
\begin{align}
g(C) &= -37.5 C/1000 - 373.41 B^{-0.6269} + 9.1256 \log(C+B) - 50.029 \\
&= -37.5 C/1000 + 9.1256 \log(C + 600) - 56.7986,
\end{align}
$$
which is the left hand side of $(1)$. Below is a plot of $g(C)$ for $0 \leq C \leq 100$:
It appears to have a root around $C = 70$, so taking $B = 600$ in $(1)$ and using Newton's method with an initial guess of $x_0 = 70$ we obtain the numerical root $C = 68.2908$. We'll call this $C_{600}$.
Similarly we find that, for $B = 601$, equation $(1)$ has a root at $C_{601} = 69.1582$.
Let's calculate our linear extrapolation formula. The line passing through the two points $(n-1,C_{n-1})$ and $(n,C_{n})$ is given by
$$
y = (C_{n}-C_{n-1})(x-n)+C_n,
$$
so taking $x = n+1$ we find that our initial guess for $C_{n+1}$ will be
$$
C_{n+1} \approx 2C_n - C_{n-1}. \tag{2}
$$
So, for example, our initial guess for using Newton's method with $B = 602$ will be
$$
C_{602} \approx 2C_{601} - C_{600} = 70.0256.
$$
Now we calculate $C_n$ iteratively for $n = 602,603,\ldots,3500$. Here is some Mathematica code which will do just that:
pts = {{600, 68.29080925572698`}, {601, 69.15821186254433`}};
lastTwo =.;
guess =.;
For[b = 602, b <= 3500, b++,
lastTwo = Take[pts, -2];
guess = 2 lastTwo[[2, 2]] - lastTwo[[1, 2]];
AppendTo[pts,
{b, c /.
FindRoot[-37.5 c/1000 - 373.41 b^(-0.6269) + 9.1256 Log[c + b] - 50.029,
{c, guess}]}
]
];
ListPlot[pts, Joined -> True]
Mathematica complains a little because of the loss of precision from using decimals instead of fractions but it still produces the desired plot:
In the Mathematica code, the FindRoot function uses Newton's method with guess as the starting point for the iteration. The variable guess is calculated in the 6th line of the code using formula $(2)$.
Best Answer
I'll try a generalization.
From $x+y-z = n $ and $x^2+y^2-z^2 = m $, $z = x+y-n$ so $z^2 =(x+y)^2-2n(x+y)+n^2 =x^2+2xy+y^2-2n(x+y)+n^2 $.
Therefore $m =x^2+y^2-(x^2+2xy+y^2-2n(x+y)+n^2) =-2xy+2n(x+y)-n^2 $ so $xy-n(x+y) =-(m+n^2)/2 $ or $(n^2-m)/2 =xy-n(x+y)+n^2 =(x-n)(y-n) $. The integer solutions thus depend on the factors of $N=(n^2-m)/2$. Note that if $n$ and $m$ have different parity, there are no integer solutions.
For each $N =ab $ either $x=n+a, y=n+b $ or $x=n-a, y=n-b $.
The largest solution is gotten by setting $a=N$ and $b = 1$ so $x = n+N$ and $y = n+1$.