Classical solutions of Neumann Laplacian I have a question concerning the Neumann Laplacian. Say, we consider a boundary value problem formulated on an interval $\Omega = [0,1]$:
\begin{equation}
-dv'' + cv = r, \quad
v'(0) = v'(1) = 0,
\end{equation}
with $c(x)\geq c_0 > 0$, $r\in L^2(0,1)$. I know, that with the help of the Lax-Milgram Lemma it is possible to show, that the problem above has a unique weak solution $v\in H^1(0,1)$. In particular, $\forall v\in H^1(0,1)$ we have 
\begin{equation}
B[v,\varphi]:= d\int_0^1 v'\varphi' dx + c\int_0^1 v\varphi dx = \int_0^1 r\varphi dx =: f[\varphi], \quad\forall \varphi\in H^1(0,1).
\end{equation}
Since the bilinear form $B[v,\varphi]$ is bounded in $H^1(0,1)$ and $H^1$-elliptic, and the functional $f[\varphi]$ is bounded as well, application of Lax-Milgram Lemma results in the existence of a unique weak solution $v\in H^1(0,1)$.
My question is: how can I make a transition to the classical solutions now? I imagine, I have to use some regularity results? 
In general, I needed to show that an operator $K:= -d\Delta + c$ is one-to-one.
Thank you for the help!
 A: Indeed, you need some regularity assumptions in order to go back from the weak
formulation using some sort of integration by parts formula. The weakest such
condition I know of is $\Delta u \in L^2$ or equivalently $\nabla u \in W^2 (
div)$ (I name this because this space appears very often in the
literature), but you could say "suppose $u \in H^2 ( \Omega)$" or even $u
\in C^2 ( \Omega)$, then prove the equivalence under this regularity. The
important thing is to be able to integrate by parts in some sense (there are
very general, abstract formulae for integration by parts).
Once you know you can do this, you take your test functions in appropriate
spaces. In your case you first take $\varphi \in C^{\infty}_0 ( \Omega)$, plug
in into the equation and use integration by parts:
$$ \int_{\Omega} r \varphi = \int_{\Omega} \nabla v \nabla \varphi +
   \int_{\Omega} v \varphi \overset{\downarrow}{=} \int_{\partial \Omega}
   \partial_n v \varphi - \int_{\Omega} \Delta v \varphi + \int_{\Omega} v
   \varphi . $$
Because of your choice for $\varphi$, the boundary term disappears. You bring
all terms together and see:
$$ \int_{\Omega} ( - \Delta v + v - r) \varphi = 0, $$
this for every $\varphi \in C^{\infty}_0 ( \Omega)$. Consequently you have the
equation
$$ - \Delta v + v = r \text{ in } L^2 ( \Omega) . $$
Now repeat the steps taking $\varphi \in C^{\infty} ( \overline{\Omega})$, and
use the equation you just found. This time the boundary term remains, but the
other disappear, from which you can deduce that the homogenenous Neumann
boundary condition holds.
Edit: Sorry, I wrote everything for dimension $N>1$ without realising that you were in $(0,1)$! You can still proceed analogously, though.
Edit: Actually, the requisite that $\Delta u \in L^2$ might be
a consequence of the first step: because we are using test functions in
$C^{\infty}_0 ( \Omega)$, we could perform the computations in the sense of
distributions, instead of using Green's formula and arrive at the equation
$$ \langle - \Delta v + v - r, \varphi \rangle_{\mathcal{D}' \times
   \mathcal{D}} = 0 $$
in the sense of distributions (as pointed out in macydanim's answer). This
would imply that $ - \Delta v = - v + r$ as distributions and (maybe under regularity assumptions on the boundary) that $\Delta
v \in L^2$, which then allows us to use Green's formula to continue.
A: You already derived a weak formulation 
\begin{align}
\int_\Omega d v' \phi ' + \int_\Omega c v \phi = \int_\Omega r \phi \quad \forall \phi \in H^1(\Omega)
\end{align}
Now suppose $v$ is sufficient smooth solution of your problem, i.e. $v\in H^1$ then it holds \begin{align}
&\int_\Omega d v' \psi ' + \int_\Omega c v \psi = \int_\Omega r \psi \quad \forall \psi \in \mathscr{D}(\Omega) \\
& \Leftrightarrow
\int_\Omega -d v'' \psi  + \int_\Omega c v \psi = \int_\Omega r \psi \quad \forall \psi \in \mathscr{D}(\Omega)
\end{align}
As $v,L,c,r \in L^2(\Omega)$, the equation $-dv''+cv=r$ is fulfilled in the distributional sense in $L^2$.
If we furthermore assume, that $v\in H^2(\Omega)$, we can use Green's formula with $u \in H^1(\Omega)$:
\begin{align}
-\int_\Omega v''u dx = \int_\Omega v'u' dx - \int_{\partial \Omega} \frac{\partial v}{\partial \nu} u d \sigma \\
\Leftrightarrow 
\int_\Omega v'u' dx  =-\int_\Omega v''u dx + \int_{\partial \Omega} \frac{\partial v}{\partial \nu} u d \sigma
\end{align}
which we can use to get 
\begin{align} 
&-d\int_\Omega v''u dx + d\int_{\partial \Omega} \frac{\partial v}{\partial \nu} u d \sigma
+\int_\Omega c v u dx= \int_\Omega r u dx  \quad &\forall u  \in H^1(\Omega)\\
 \Leftrightarrow &
-d\int_\Omega v''u dx + d\int_{\partial \Omega} \frac{\partial v}{\partial \nu} u d \sigma
+\int_\Omega c v u dx- \int_\Omega r u dx =0  \quad &\forall u  \in H^1(\Omega)\\
 \Leftrightarrow &\underbrace{\int_\Omega (-dv'' +cv -r) u dx}_{=0\text{, as v is a weak solution}} +d\int_{\partial \Omega} \frac{\partial v}{\partial \nu} u d \sigma=0  \quad &\forall u  \in H^1(\Omega)\\
 \Leftrightarrow & \int_{\partial \Omega} \frac{\partial v}{\partial \nu} u d \sigma=0 \quad &\forall u  \in H^1(\Omega)
\end{align}
And since the trace space on $\Omega$ is dense in $L^2(\Omega)$:
 \begin{align}
\int_{\partial \Omega} \frac{\partial v}{\partial \nu} w d \sigma=0 \quad &\forall w  \in L^2(\Omega)
\end{align}
So also the Neumann boundary condition are fulfilled in a $L^2$ sense.
Reference is chapter 4.2.
