Thanks for your help everyone. I think I got it. Kindly check my proof for errors if any:
Proof:
$\forall x \in [0,1]$, fix $\epsilon > 0$. Now because at every point, both left and right hand limits exist(except at $0$ and $1$, in which case only one of the following will be there),
$$\exists \delta_x,\eta_x > 0$$
$$x_0 \in (x-\delta_x,x) \Rightarrow |f(x_0) - f(x^-)| < \epsilon$$
$$x_0 \in (x,x+\eta_x) \Rightarrow |f(x_0) - f(x)| < \epsilon$$
where $f(x^-)$ is the left hand limit of the function at x and since we have right continuity, $f(x^+) = f(x)$.
Let $O_x = (x-\delta_x,x)$ and $V_x = (x,x+\eta_x)$. Let $W_x = O_x\cup V_x \cup\{x\}$, which is still an open set. Then
$$\bigcup_{x \in [0,1]} W_x$$
is an open cover for $[0,1]$. Since $[0,1]$ is compact, it has a finite subcover. Hence
$$ [0,1] \subseteq \bigcup_{k=1}^n W_{x_k}$$
For $k=1...n$, if $x \in O_{x_k}$, then $|f(x)-f(x^-_k)| < \epsilon \Rightarrow |f(x)| < |f(x^-_k)| + \epsilon$
If $M_1 = \max_{k=1...n} |f(x_k^-)|$, then $|f(x)| < M_1 + \epsilon$. (Only if $x$ belongs to any of the "$O$" balls)
Similarly, by defining $M_2 = \max_{k=n+1...m} |f(x_k)|$, and following similar steps, we get $|f(x)| < M_2 + \epsilon$. (For x belonging to "$V$" balls).
Finally, if $x \in \{x_1,\cdots, x_n\}$, then $|f(x)| \le M_3 = \max_{1\le i\le n} |f(x_i)|$.
So the bound for $x \in [0,1]$ is,
$$|f(x)| < \max(M_1,M_2,M_3) + \epsilon$$
Thus f is bounded. $\blacksquare$
I am interested in other elegant proofs if there are any. Thanks again.
Update: Thanks to $\epsilon-\delta$ for the comment. I have fixed the proof, I think. Please check it.
If $f$ is unbounded near $a$, it can clearly not have bounded variation, since
$$V(f,[a,b]) \geqslant \lvert f(x)-f(a)\rvert$$
for all $x \in (a,b]$.
If $-\infty < m = \liminf\limits_{x \searrow a} f(x) < \limsup\limits_{x\searrow a} f(x) = M < +\infty$, let $\varepsilon = (M-m)/3$. Then we can find sequences $x_n \searrow a$ and $y_n \searrow a$ with $x_n > y_n > x_{n+1}$ and $f(y_n) < m+\varepsilon$, $f(x_n) > M-\varepsilon$ for all $n$, so
$$V(f,[a,b]) \geqslant \sum_{n=1}^N \lvert f(x_n) - f(y_n)\rvert \geqslant N\cdot \varepsilon.$$
Letting $N\to \infty$, we see that then $V(f,[a,b]) = +\infty$.
Best Answer
the function is discontinuous for all rational numbers, which results in problems with interpreting and working with the function. However, the limit there is zero. That means those points are removable discontinuities, which makes the popcorn function (weirdly) a continuous function, even though it isn’t intuitively so.