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.