Applying Schwinger-Dyson equations within the LSZ formula

My problem will be formulated in terms of $$\phi^3$$ theory, and I would appreciate answers within the framework of $$\phi^3$$ or another scalar field theory. This question is to help me understand what people mean by "contact terms can be ignored when calculating scattering amplitudes." I have often heard this in the context of deriving $$k_\mu\mathcal M^\mu=0$$ in the context of QED. There have been a healthy number of questions in that QED context about justifying that relation and whether contact terms vanish in the LSZ formula or instead survive but fail to contribute connected diagrams. I have constructed a toy example here that illustrates my personal confusion about contact terms.

The Lagrangian density is $$\mathcal L= -\frac{1}{2}Z_{\phi}\partial^{\mu}\phi\partial_{\mu}\phi -\frac{1}{2}Z_mm^2\phi^2+\frac{1}{6}Z_gg\phi^3+Y\phi$$

I want to look at the $$\phi \phi \rightarrow \phi\phi$$ scattering amplitude at the tree level of $$g^2$$. We find this amplitude by the LSZ formula:

$$\langle f| i\rangle = i^4\int dx_1 dx_2 dx_1' dx_2' e^{i(k_1x_1+k_2x_2-k_1'x_1'-k_2'x_2')}(-\partial^2_1+m^2)(-\partial^2_2+m^2)(-\partial^2_3+m^2)(-\partial^2_4+m^2)\langle0|T\phi(x_1)\phi(x_2)\phi(x_1')\phi(x_2')|0\rangle$$

Now, before rewriting that correlation function in terms of the free-vacuum and using Wick's theorem/writing the correlation function as functional derivatives of some generating function, let's instead first apply the Schwinger-Dyson equations.

First, note that for the Lagrangian above, the Schwinger-Dyson equations imply that for some correlation function with n fields, $$(-Z_\phi \partial_{x_1}^2+Z_m m^2-Y)i\langle0|T\phi(x_1)...\phi(x_n)|0\rangle = \frac{1}{2}g\langle0|T\phi(x_1)^2...\phi(x_n)|0\rangle+\sum_{j=2}^n\delta(x-x_j)\langle0|T\phi(x_2)...\phi(x_{j-1})\phi(x_{j+1})...\phi(x_n)|0\rangle.$$

The terms in the sum including delta-functions are typically called contact terms.

Now I imagine a small devil on my shoulder telling me that contact terms may be neglected as they fail to produce the right singularities. This appears to me to be the accepted answer here for example. Heeding his advice, I apply the Schwinger-Dyson equations four times, one for each differential operator in the LSZ formula. At tree-level, I'll take the $$Z$$ factors as 1 and $$Y$$ as 0.

$$\langle f| i\rangle = g^4\int dx_1 dx_2 dx_1' dx_2' e^{i(k_1x_1+k_2x_2-k_1'x_1'-k_2'x_2')}\frac{1}{16}\langle0|T\phi(x_1)^2\phi(x_2)^2\phi(x_1')^2\phi(x_2')^2|0\rangle + \text{graveyard of neglected contact terms}$$

The devil led me astray, since the scattering amplitude should go as $$g^2$$, not $$g^4$$! This now leads me at last to my question, in two parts:

1. What went wrong in the example above? Namely, was it that I could not neglect contact terms after all? Or was there another mistake which when rectified would allow me to drop all contact terms?

2. Could you show me how to get the correct tree-level $$g^2$$-proportional scattering using the Schwinger-Dyson equations at least once within the LSZ formula?