In your proposed solution, all of your rectangles have the same width. The key is to let the width of the rectangles vary which will make the summation significantly easier.
Let $x_i = \frac{i^2}{n^2}$ for all $n \leq i \leq 3n$. Notice that $f(x_n) = 1$ and $f(x_{3n}) = 9$. Then $$\Delta x_i = x_{i+1} - x_i= \frac{(i+1)^2}{n^2} - \frac{(i)^2}{n^2} = \frac{2i + 1}{n^2}.$$
Then the (left) Riemann integral becomes $$\lim_{n \to \infty} \sum_{i = n}^{3n -1} f(x_i) \Delta x_i = \lim_{n \to \infty} \sum_{i = n}^{3n-1} \frac{i}{n} \frac{2i + 1}{n^2},$$ which you can calculate using the standard formula $$\sum_{i = 1}^n i^2 = \frac{n(n+1)(2n+1)}{6}.$$
I will remark that of all the rectangle widths, the largest is $\Delta x_{(3n-1)}$, where we have $$\Delta x_{(3n-1)}=\frac{(3n)^2 - (3n-1)^2}{n^2} \to 0,$$ as $n \to \infty$. Thus all the widths of the rectangles are becoming arbitrarily small.
Let
$$
x_n=\sum_{k=0}^nf\left(\frac{k}{n}\right)-n\int_0^1f(x)dx
$$
We will use the following result:
Lemma If $g:[0,1]\to\mathbb{R}$ is a continuously differentiable function. Then
$$
\frac{g(0)+g(1)}{2}-\int_0^1g(x)dx= \int_{0}^1\left(x-\frac{1}{2}\right)g'(x)dx.
$$
Indeed, this is just integration by parts:
$$\eqalign{
\int_{0}^1\left(x-\frac{1}{2}\right)g'(x)dx
&=\left.\left(x-\frac{1}{2}\right)g(x)\right]_{x=0}^{x=1}
-\int_0^1g(x)dx\cr
&=\frac{g(1)+g(0)}{2}-\int_0^1g(x)dx
}$$
Now applying this to the functions
$x\mapsto f\left(\frac{k+x}{n}\right)$ for $k=0,1,\ldots,n-1$ and adding the resulting inequalities we obtain
$$
x_n-\frac{f(0) +f(1)}{2} = \int_0^1\left(x-\frac{1}{2}\right)H_n(x)dx\tag{1}
$$
where,
$$
H_n(x)=\frac{1}{n}\sum_{k=0}^{n-1}f'\left(\frac{k+x}{n}\right)
$$
Clearly for every $x$, $H_n(x)$ is a Riemann sum of the function continuous $f'$, hence
$$
\forall\,x\in[0,1],\quad\lim_{n\to\infty}H_n(x)=\int_0^1f'(t)dt
$$
Moreover, $| H_n(x)|\leq\sup_{[0,1]}|f'|$. So,
taking the limit in $(1)$ and
applying the Dominated Convergence Theorem, we obtain
$$
\lim_{n\to\infty}\left(x_n-\frac{f(0) +f(1)}{2}\right)=
\left(\int_0^1f'(t)dt\right)\int_0^1\left(x-\frac{1}{2}\right)dx=0.
$$
This proves that
$$
\lim_{n\to\infty}x_n=\frac{f(0) +f(1)}{2}
$$
And consequently
$$
\lim_{n\to\infty}\left(\sum_{k=\color{red}{1}}^nf\left(\frac{k}{n}\right)-n\int_0^1f(x)dx\right)=\frac{f(1)-f(0)}{2}
$$
Best Answer
Let $x_k=3+\dfrac{9k}n$ and $\Delta x=\dfrac9n$, so
$$\sum_{k=1}^n f(x_k)\Delta x=\sum_{k=1}^n \left(3+\dfrac{9k}n\right)^2\dfrac9n=\sum_{k=1}^n \left(9+\dfrac{54k}n+\dfrac{81k^2}{n^2}\right)\dfrac9n$$
$$=\left(9n+\dfrac{54}n\dfrac{n(n+1)}2+\dfrac{81}{n^2}\dfrac{n(n+1)(2n+1)}6\right)\dfrac9n,\\$$
which approaches $(9+27+27)9=567$ as $n$ approaches $\infty$.
Your answer also approaches $567$ as $n\to\infty,$ but you had $-$ where I have $+$.