Integrating by parts, using a cutoff function if necessary, we have
$$ \int_{B(r)} |\nabla u|^2\ll \int_{B(r)} a_{ij}\partial_iu\partial_ju=\int_{B(R)} uLu+O(u|\nabla u|)\ll \int_{B(R)} |u\nabla u|. $$
Using Cauchy Schwarz, you can get
$$ \int_{B(r)} |\nabla u|^2 \ll \int_{B(R)} |u|^2. $$
(If $u$ does not have $W^{1,2}$ regularity to start with, you can approximate the left hand side by difference quotients.)
Now do the same for $\nabla u$, $\nabla^2 u$, etc. You will get
$$ \|u\|_{W^{k,2}(B(1/2))}\ll \|u\|_{L^2(B(1))} $$
for each $k$. Choose a large enough $k$ so that you have
$$ \|u\|_{W^{2,\infty}} \ll \|u\|_{W^{k,2}} $$
and you're done.
NOTE: The following solution that I originally posted is wrong.
This is an expansion of Yang's suggestion. Moser iteration (combined with some Maximum-principle style arguments, see Chapter 8 of Gilbarg-Trudinger) would give you $C^\alpha$ bound for $u$. Then you fix the coefficients locally and rewrite the PDE as
$$\partial_tu+\Delta u=(\Delta-L)u$$
and treat the right hand side as an inhomogeneous term (which has been shown to lie in $C^\alpha_{loc}$). Then Schauder theory for the Poisson equation (Chapter 4 of Gilbarg-Trudinger) would give you local $C^{2,\alpha}$ bound for $u$, which is more than what you have asked.