is not matching
Because it does not "look" the same, does not mean it is not equivalent or that it is wrong.
Mathematica solution is correct, it is just written differently. One way to verify, is to simply pluging some values for s and F[x] and see if Mathematica solution matches the Green function solution by hand. So below I solved the problem using Green Function method, and compared the solution for say s=2 and f(x)=x and got the same plot.
Mathematica seems to just expanded the integrals out, that come from the Green function integrals.
ClearAll[f,s,x,y]
mmaSol=y[x]/.First@DSolve[{y''[x]-s*y[x]==-f[x],y'[0]==0,y'[1]==0},y[x],x]
Plot[Evaluate[mmaSol/.{s->2,f[z_]->z}],{x,0,1}]
Now I found the Green function solution, and plotted it using same force and s and got same exact solution
solve
\begin{align*}
y^{\prime\prime}\left( x\right) -sy\left( x\right) & =-f\left(
x\right) \\
y^{\prime}\left( 0\right) & =0\\
y^{\prime}\left( 1\right) & =0
\end{align*}
The first step is to determine $y_{1},y_{2}$. These are the two fundamental
solutions of $y^{\prime\prime}-sy=0$ . To simplify the derivation, let
$y_{1}\left( x\right) $ be the solution that satisfies the boundary
conditions at the left end of domain and $y_{2}\left( x\right) $ satisfies
the boundary condition on the right end.
The solution to $y^{\prime\prime}\left( x\right) -sy\left( x\right) =0$ is
$y_{h}=c_{1}e^{\sqrt{s}x}+c_{2}e^{-\sqrt{s}x}$. Hence $y_{h}^{\prime}
=c_{1}\sqrt{s}e^{\sqrt{s}x}-c_{2}\sqrt{s}e^{-\sqrt{s}x}$. Applying
BC $y^{\prime}\left( 0\right) =0$ gives
\begin{align*}
0 & =c_{1}\sqrt{s}-c_{2}\sqrt{s}\\
c_{1} & =c_{2}
\end{align*}
Therefore
\begin{align*}
y_{1}\left( x\right) =e^{\sqrt{s}x}+e^{-\sqrt{s}x}
\end{align*}
Applying BC $y^{\prime}\left( 1\right) =0$ gives
\begin{align*}
0 & =c_{1}\sqrt{s}e^{\sqrt{s}}-c_{2}\sqrt{s}e^{-\sqrt{s}}\\
c_{2} & =c_{1}\frac{\sqrt{s}e^{\sqrt{s}}}{\sqrt{s}e^{-\sqrt{s}}}\\
& =c_{1}e^{2\sqrt{s}}
\end{align*}
Therefore
\begin{align*}
y_{2}\left( x\right) =e^{\sqrt{s}x}+e^{2\sqrt{s}}e^{-\sqrt{s}x}
\end{align*}
The Wronskian is
\begin{align*}
W & =
\begin{vmatrix}
y_{1} & y_{2}\\
y_{1}^{\prime} & y_{2}^{\prime}
\end{vmatrix}
=
\begin{vmatrix}
e^{\sqrt{s}x}+e^{-\sqrt{s}x} & e^{\sqrt{s}x}+e^{2\sqrt{s}}e^{-\sqrt{s}x}\\
\sqrt{s}e^{\sqrt{s}x}-\sqrt{s}e^{-\sqrt{s}x} & \sqrt{s}e^{\sqrt{s}x}
-e^{2\sqrt{s}}\sqrt{s}e^{-\sqrt{s}x}
\end{vmatrix}
\\
& =\left( e^{\sqrt{s}x}+e^{-\sqrt{s}x}\right) \left( \sqrt{s}e^{\sqrt{s}
x}-e^{2\sqrt{s}}\sqrt{s}e^{-\sqrt{s}x}\right) -\left( e^{\sqrt{s}
x}+e^{2\sqrt{s}}e^{-\sqrt{s}x}\right) \left( \sqrt{s}e^{\sqrt{s}x}-\sqrt
{s}e^{-\sqrt{s}x}\right) \\
& =2\sqrt{s}\left( 1-e^{2\sqrt{s}}\right)
\end{align*}
For $s\neq0$ the above is non-zero. The Green function is
\begin{align*}
G\left( x,a\right) & =\frac{-1}{W}\left\{
\begin{array}
[c]{ccc}
y_{1}\left( a\right) y_{2}\left( x\right) & & 0\leq a\leq x\\
y_{1}\left( x\right) y_{2}\left( a\right) & & x\leq a\leq1
\end{array}
\right. \\
& =\frac{1}{2\sqrt{s}\left( e^{2\sqrt{s}}-1\right) }\left\{
\begin{array}
[c]{ccc}
\left( e^{\sqrt{s}a}+e^{-\sqrt{s}a}\right) \left( e^{\sqrt{s}x}
+e^{2\sqrt{s}}e^{-\sqrt{s}x}\right) & & 0\leq a\leq x\\
\left( e^{\sqrt{s}x}+e^{-\sqrt{s}x}\right) \left( e^{\sqrt{s}a}
+e^{2\sqrt{s}}e^{-\sqrt{s}a}\right) & & x\leq a\leq1
\end{array}
\right.
\end{align*}
Hence the solution is
\begin{align*}
y\left( x\right) =\int_{0}^{x}G\left( x,a\right) F\left( a\right)
da+\int_{x}^{1}G\left( x,a\right) F\left( a\right) da
\end{align*}
Or
\begin{align*}
y\left( x\right) & =\frac{1}{2\sqrt{s}\left( e^{2\sqrt{s}}-1\right)
}\int_{0}^{x}\left( e^{\sqrt{s}a}+e^{-\sqrt{s}a}\right) \left( e^{\sqrt
{s}x}+e^{2\sqrt{s}}e^{-\sqrt{s}x}\right) F\left( a\right) da\\
& +\frac{1}{2\sqrt{s}\left( e^{2\sqrt{s}}-1\right) }\int_{x}^{1}\left(
e^{\sqrt{s}x}+e^{-\sqrt{s}x}\right) \left( e^{\sqrt{s}a}+e^{2\sqrt{s}
}e^{-\sqrt{s}a}\right) F\left( a\right) da
\end{align*}
Mathematica result is an expansion of the above integral. Mathematica also using d[K1] and d[K2] for the integrals, while in the hand solution da is used for both. But it is the same.
Plotting the above
ClearAll[x,y,arg,y1,y2,s,f];
y1[arg_]:=Exp[Sqrt[s] arg]+Exp[-Sqrt[s] arg];
y2[arg_]:=Exp[Sqrt[s]arg]+Exp[2 Sqrt[s]] Exp[-Sqrt[s] arg];
wronskian=Det[{{y1[x],y2[x]},{D[y1[x],x],D[y2[x],x]}}]//Simplify
mySol=-1/wronskian Integrate[y1[a] y2[x] f[a],{a,0,x}]
-1/wronskian Integrate[y1[x] y2[a] f[a],{a,x,1}]
Plot[Evaluate[mySol/.{s->2,f[z_]->z}],{x,0,1}]
You can try different s values and different f. You'll get same result as Mathematica's. Sometimes it is hard for CAS to give same looking result as hand solution.
