Suppose that $f:[0,\infty)\to[-1,1]$ is a continuous function at zero and at all but finite number of points such that $f(0)=1$ and $\int_0^\infty f^2(x)dx<\infty$. Let $\{a_n\}_{n\ge1}$ be a sequence of positive integers such that $a_n\to\infty$ as $n\to\infty$. The right Riemann sum $$ \frac1{a_n}\sum_{j=1}^{a_n}f^2(j/a_n) $$ converges to the integral $\int_0^1f^2(x)dx$ as $n\to\infty$. Also, for a positive integer $m$, we have that $$ \frac1{a_n}\sum_{j=1}^{ma_n}f^2(j/a_n)\to\int_0^mf^2(x)dx\quad\text{as}\quad n\to\infty. $$ Let $\{b_n\}_{n\ge1}$ be a sequence of positive integers such that $b_n/a_n\to\infty$ as $n\to\infty$. Intuitively, it seems that $$ \frac1{a_n}\sum_{j=1}^{b_n}f^2(j/a_n)\to\int_0^\infty f^2(x)dx\quad\text{as}\quad n\to\infty. $$ Is it true? How can I show this?
References are also welcome. Any help is much appreciated!

I'm not sure why you square $f.$ Why not just assume $f:[0,\infty)\to [0,1]$ continuously, and $\int_0^\infty f(x)\,dx < \infty?$
A counterexample given by Alex Francisco shows the answer to your question is no. I'm going to describe another answer that uses the same idea, but may be simpler.
Over each $n=2,3,\dots$ consider a thin triangle centered over $n,$ of base $1/n^2$ and height $1.$ Let $f$ be the tops of these triangles, with $f=0$ everywhere else. Then $f$ is continuous on $[0,\infty),$ and
$$\int_0^\infty f(x)\, dx = \sum_{n=2}^{\infty}\left (\frac{1}{2}\cdot\frac{1}{n^2}\cdot 1\right ) <\infty.$$
However, thinking of $a_n=n, b_n = n^3,$ we have
$$\frac{1}{n}\sum_{j=1}^{n^3}f(j/n) \ge \frac{1}{n}(f(2) + f(3) + \cdots +f(n^2)) = \frac{1}{n}(n^2-1) \to \infty.$$