Solution (source code)

= Solution

Continuity at $x=\xi$ is automatic in the proposed expression. Integrating the differential equation through $x=\xi$ gives the <derivative> jump
$$
G_x(\xi^+,\xi)-G_x(\xi^-,\xi)=1.
$$
Thus
$$
c(\xi)\left[y_1(\xi)y_2'(\xi)-y_1'(\xi)y_2(\xi)\right]=1.
$$
In terms of the <Wronskian>
$$
W(x)=y_1(x)y_2'(x)-y_1'(x)y_2(x),
$$
one has
$$
\boxed{c(\xi)=\frac1{W(\xi)}}.
$$
The <boundary conditions> hold because $y_1'(0)=0$ and $y_2'(1)=0$.

The <Abel identity> gives $W'=-\alpha W$. If $\alpha=0$, the Wronskian and hence $c$ are constant. The two branches of the displayed formula are then interchanged by $x\leftrightarrow\xi$, proving
$$
\boxed{G(x,\xi)=G(\xi,x)}.
$$
This is the symmetry of the <Neumann Green function for a second-order ordinary differential equation> in the self-adjoint case.

For $G''-G=\delta(x-\xi)$, choose
$$
y_1(x)=\cosh x,
\qquad
y_2(x)=\cosh(1-x).
$$
Their Wronskian is
$$
W=-\cosh x\sinh(1-x)-\sinh x\cosh(1-x)=-\sinh1.
$$
Writing $x_<=\min(x,\xi)$ and $x_>=\max(x,\xi)$ gives
$$
\boxed{G(x,\xi)=-\frac{\cosh x_<\cosh(1-x_>)}{\sinh1}}.
$$
The solution of the inhomogeneous problem is $y(x)=\int_0^1G(x,\xi)\xi\,d\xi$. Equivalently, solving directly gives
$$
y(x)=A\cosh x+B\sinh x-x.
$$
The two Neumann conditions yield $B=1$ and
$$
A=\frac{1-\cosh1}{\sinh1}=-\tanh\frac12.
$$
Therefore
$$
\boxed{y(x)=\sinh x-\tanh\left(\frac12\right)\cosh x-x}.
$$