I am currently trying to prove that for all $\epsilon >0$ and $f \in \mathcal{L}^1(\mathbb R^n,\mathbb K)$ where $\mathbb K$ is either $\mathbb R$ or $\mathbb C$ that there is a continuous function $g$ such that $\int |f-g| d\mu <\epsilon$ where $\mu$ is the Lebesgue measure.
I already proved that for all $\epsilon >0$ there is a simple function $\psi=\sum_{j=1}^{m}a_j1_{E_j}$ where $E_j$ is a finite union of open balls and $1_{E_j}$ is the characteristic function on $E_j$ such that $\int |f-\psi| d\mu < \epsilon$.
Hence the theorem is easily proved once I show that each $1_{E_j}$ can be approximated by continuous functions, in the sense that for all $\epsilon >0$ there is a continuous function $g$ such that $\int |1_{E_j}-g| d\mu <\epsilon$.
For $n=1$, this is certainly an easy task, as I can construct a sequence $(g_n)$ with $g_n|_{E_j}=1$ and $g_n$ decreasing to $0$ linearly outside of $E_j$ such that $\lim g_n = 1_{E_j}$. However, I don't see how this can be proved for $n>1$ and it doesn't seem like there is a generalization of the approach for the case of $n=1$.
One way to do that in $\mathbb{R}^n$ is to use the distance to $E_j$. We define for $x \in \mathbb{R}^n$ and $A \subset \mathbb{R}^n$ :
$$d(x,A):= \inf \{ \| x - y\| \mid y \in A\}.$$
The function $x \mapsto d(x,A)$ is continuous (one can see that it is actually 1-Lipschitz) and is equal to $0$ on $A$ (and on its topological closure). Now, let $E'_j$ be the obtained from $E_j$ by doing the union of the same open balls but with a radius increased by $\varepsilon$ for each ball and let $F_j = \mathbb{R}^n \setminus E_j$. You can imagine $E_j'$ as a slightly bigger $E_j$, and $F_j$ as the complement of this : you then have $E_j$, a small "annulus" around $E_j$, and everything else is in $F_j$.
We now define :
$$g :x \mapsto \dfrac{d(x,F_j)}{d(x,F_j)+d(x,E_j)}$$
This function goes into $[0,1]$ ; it is equal to $1$ on $E_j$, gets smaller as $x$ gets further from $E_j$ and closer to $F_j$, and is equal to $0$ on $F_j$. The integral of $|f_j - g|$ on $\mathbb{R}^n$ is therefore $$\begin{array}{lll} \int_{\mathbb{R}^n} |f_j - g|= \int_{F_j} |f_j - g|+ \int_{E_j} |f_j - g|+ \int_{\mathbb{R}^n \setminus (E_j \cup F_j)} |f_j - g| \\ = \int_{F_j} |0 - 0|+ \int_{E_j} |1-1|+ \int_{\mathbb{R}^n \setminus (E_j \cup F_j)} |0 - g| \\ = \int_{\mathbb{R}^n \setminus (E_j \cup F_j)} |g| \\ \end{array}$$ As $g$ has its values in $[0,1]$, we thus have $$\int_{\mathbb{R}^n} |f_j - g| \leqslant \int_{\mathbb{R}^n \setminus (E_j \cup F_j)} 1$$ Given the construction of $F_j$ from $E_j$, the volume of $\mathbb{R}^n \setminus (E_j \cup F_j)$ is less than $C\times \varepsilon$, where $C$ is a constant. Thus $$\int_{\mathbb{R}^n} |f_j - g| \leqslant C \times \varepsilon$$ and you can conclude by taking $\varepsilon$ small enough.