Let's follows Evans's PDE book notation. We say $G(x,y)$ is the green function build on $\Omega\subset \mathbb R^N$ where $\Omega$ open bounded nice boundary.
Then if $u\in C^2(\bar\Omega)$, we could represent $u$ by $$ u(x)=\int_\Omega G(x,y)\Delta u(y)\,dy+\int_{\partial\Omega}\partial_\nu G(x,y)u(y)\,d\sigma(y) \tag 1$$
Ok, for $u\in C^2(\bar \Omega)$, proving $(1)$ is an easy task. However, my professor states that this formula is also true if we only have $u\in C^2(\Omega)\cap C^0(\bar\Omega)$ and we could prove this by using approximation...
I tried a while but I still got stuck on this approximation problem.
Honestly, I don't even have a good idea how to approx a function $u\in C^2(\Omega)\cap C^0(\bar\Omega)$ by $(u_n)\subset C^2(\bar\Omega)$.
I think the approximation should hold in $L^\infty$ sense so I was trying to using mollification... But I don't really see how...Any help is really welcome!
Here @Tomas suggest we go approx with respect to domain. Here is what I tried (didn't get too much however...). Define $$\Omega_\epsilon:=\{x\in\Omega,\,\,\text{dist}(x,\partial\Omega)>\epsilon\} $$ Now fix arbitrary $x\in \Omega$ and for $\epsilon>0$ small enough we have $x\in\Omega_\epsilon$. Also we have $u\in C^2(\bar \Omega_\epsilon)$. Apply green identity, we have $$ u(x)=\int_{\Omega_\epsilon}\Delta u(y)G(x,y)dy+\int_{\partial\Omega_\epsilon}\partial_\nu G(x,y)u(y)\,d\sigma - \int_{\partial\Omega_\epsilon} G(x,y)\partial_\nu u(y)\,d\sigma:=I+II+III$$
The Problem is how we take $\epsilon\to 0$. We have no problem with $II$, but $I$ and $III$ will be the problem.... In order to much to limit by using LDCT, we need to know, for $I$ $$ \int_{\Omega}\Delta u(y)G(x,y)dy $$ is well define. And hence we could use LDCT to justify $I$. Of course if $u\in C^2(\bar\Omega)$ we are good to go; if not, we have to consider the behavior of $\Delta u$ around $\partial \Omega$. As $G(x,y)=0$ over $\partial\Omega$ we know that $G(x,y)\to 0$ as $y\to \partial\Omega$. Hence $G(x,y)\to 0$ uniformly. Now the problem is $G(x,y)$ compete with $\Delta u$ around the boundary...but could $\Delta u$ just blow up at $\partial\Omega$ and hence the integration blows up?
The part $III$ actually have the same problem as part $I$...
I think the best you can do is the following: There exists $p_n$ such that for every $u\in W^{2,p}(\Omega)$ with $p>p_n$ Green's representation formula is valid for $u$. To see this note that for $2p>N$ $u$ is continuous in $\bar{\Omega}$ by Sobolev embedding, so that the LHS and the boundary term make sense. Moreover, since $G(x,\cdot)\in L^q(\Omega)$, uniformly in $x$, for $1<q<N/(N-2)$ we get that for $p>N/2$ the whole formula makes sense. The fact that $C^\infty(\bar{\Omega})$ is dense in $W^{2,p}(\Omega)$ follows since $\Omega$ is smooth.
To see that this is optimal consider $y_0\in \partial \Omega$ and consider $u(x)=|x-y_0|^2\sin(|x-y_0|^{2a})$. For values $a<<-2N$, we get that $u\in C^2(\Omega) \cap C(\bar{\Omega})$, but the first integral term in the representation formula doesn't make sense, since $\Delta u$ blows up too fast near $y_0$. If by $u\in C^2(\Omega)$ you mean that $D^2 u$ is bounded, then interpolation inequalities give that $u\in W^{2,p}(\Omega)$ for every finite $p$, and thus we're in the situation above.
Of course all of this handwaves the approximation, citing a more general result for Sobolev spaces. I don't know if you can prove a similar result avoiding this, and constructing a direct approximating sequence...