I was trying to write the solution of an inhomogeneous differential equation $(\partial^2_x+m^2)\phi(x)=\rho(x)$ using the Green function: \begin{equation} (\partial_x^2+m^2)G(x,y)=\delta(x-y). \end{equation} Then my solution is \begin{equation} \phi(x)=\phi_0(x)+\int dy~ \rho(y)G(x,y), \end{equation} where $\phi_0(x)$ is the solution for the homogeneous case (i.e. when $\rho=0$).
Then I asked myself if it was possible to write everything in form of an integral with the Green function multiplying the whole argument. Something like $\phi=\int dyG(x,y)(...)$.
As first step I used a $\delta$-function on $\phi_0$, then from the definition of the Green function I translated the $\delta$-function in terms of $G(x,y)$. As last step I integrated by part (assuming that $\phi(x)$ is a local field and it goes to zero at infinity) \begin{align} \phi(x)&=\int dy~\{((\partial_y^2+m^2)G(x,y))\phi_0(y)+\rho(y)\} \\ &=\int dy ~G(x,y)\{(\partial_y^2+m^2)\phi_0(y)+\rho(y)\}. \end{align} Problem: With this last relation the first term disappears, since $\phi_0$ is a solution of the homogeneous case. Hence \begin{equation} \phi(x)=\int dy~G(x,y)\rho(y), \end{equation} which is wrong in general, since the term satisfying the homogeneous case is now missing.
Question: With $D_x:=(\partial_x^2+m^2)$. What's wrong in using the Green's function defined as $ D_x G(x,y) =\delta (x-y) $, to write: $$ D_x \phi_0(x)=0,\qquad \phi_0(x)=\int dy \phi_0(y)\delta (x-y)= \int dy \phi_0(y)D_yG (x,y)= \int dy G (x,y)D_y\phi_0(y)=0?$$
Any suggestions?
The problem you're having is that in the case where the boundary conditions are homogenous (i.e., $\phi$ vanishes at the boundaries), the homogeneous solution is identically zero. Which you've proven, so yay! But let's see what happens when you try to use your technique to solve a differential equation with different boundary conditions, where the homogeneous solution shouldn't vanish. What we're going to find is that in this case, the boundary terms you get when you integrate by parts no longer vanish, and so the homogeneous function is no longer zero.
In particular, let's try to solve this equation on the domain $[x_1,x_2]$ when $\phi = a \neq 0$ on the boundary. (Let's also, for the sake of argument, work in 1-D. I'm not sure if that's what you intended, but it makes the argument easiest.) In this case, the solution to your equation can be written as $$ \phi(x) = \phi_0 (x) + \int G(x,y) \rho(y) \, dy $$ where $G(x,y)$ satisfies the homogeneous boundary conditions as before: $$ (\partial_x^2 + m^2) G(x,y) = \delta(x-y), \quad G(x_1,y) = G(x_2,y) = 0 $$ and $\phi_{in}$ is the solution to the homogenous equation with inhomogeneous boundary conditions: $$ (\partial_x^2 + m^2) \phi_0 = 0, \quad \phi_0(x_1) = \phi_0(x_2) = a. $$ It is not too hard to show that this expression for $\phi(x)$ satisfies the desired differential equation with the boundary conditions $\phi = a$ at $x_1$ and $x_2$.
Now let's try to use your logic to find the homogeneous solution $\phi_0$. If we do this, we get \begin{align} \phi_0(x) &= \int \phi_0(y) \delta(x - y) \, dy \\ &= \int \phi_0(y) (\partial_y^2 + m^2) G(x,y) \, dy \\ &= \int G(x,y) (\partial_y^2 + m^2) \phi_0(y) \, dy + \left[ \phi_0(y) \partial_y G(x,y) \right]_{y=x_1}^{x_2} - \left[ (\partial_y \phi_0) G(x,y) \right]_{y=x_1}^{x_2} \end{align} Now, the integral vanishes by the logic you outlined in your question, and the second boundary term vanishes by the boundary conditions we've imposed on $G(x,y)$ (and the fact that it's symmetric.) But the first term does not in fact vanish, and in fact is equal to $$ \phi_0(x) = a \left[ \frac{\partial G(x,y)}{\partial y} \right]_{y=x_1}^{x_2}. $$ We conclude that the homogeneous term does not vanish in this case—unless $a = 0$, in which case it does vanish, just as you found. The point is, though, that you can't in general discard the boundary conditions if you're trying to solve the problem with inhomogeneous boundary conditions; and they will give rise to a non-zero homogeneous solution.
Just to convince you that I'm not pulling a fast one on you: let's look at the case $-x_1 = x_2 = 1$. In this case, the Green's function with homogenous boundary conditions is $$ G(x,y) = \frac{1}{m \sin(2 m)}\begin{cases} \sin(m(y-1)) \sin (m(x+1)) & x< y \\ \sin(m(y+1)) \sin (m(x-1)) & x > y \end{cases} $$ Its derivative with respect to $y$ is $$ \frac{\partial G(x,y)}{\partial y} = \frac{1}{\sin(2m)}\begin{cases} \cos(m(y-1)) \sin (m(x+1)) & x< y \\ \cos(m(y+1)) \sin(m(x-1)) & x > y \end{cases} $$ By our above derivation, we need to evaluate this at $y = \pm 1$. When $y = 1$, we have $x < y$, and so we're in the first case above; when $y = -1$, we're in the second. In both cases, the $y$-dependent terms are $\cos(0) = 1$. All told, then, $$ \phi_0 = a \frac{\sin (m(x+1)) - \sin(m (x-1))}{\sin (2m)} = a \frac{2 \sin(m)}{\sin(2m)} \cos (mx) = a \frac{\cos{mx}}{\cos m}. $$ Thus, the homogeneous solution is a cosine curve that goes to $a$ when $x = \pm 1$, as we would have expected.