this is my first question on Stack-Exchange. I am not an expert in numerical integration and I am having issues of convergence with some numerical integration. I am trying to solve numerically a nonlinear differential equation of which I know an exact solution. I want to compare my numerical results with the exact one (to check the correctness of the code). Unfortunately, after 2 iterations the numerical result starts departing from the exact behavior.

This is the differential equation:

y”(t)=(y(t)^2-1)y(t)

which I wrote in an integral form (I need to solve it in the integral form for other reasons) in terms of the Green function

G(t-x)=(1/2) (t-x) sign(t-x).

An analytic solution is given by

y(x)=Tanh(x/sqrt{2})

The Mathematica code I am using is the following:

```
t1 = -10;
t2 = 10;
n = 100;
m = 5;
t = Range(t1, t2, (t2 - t1)/(n));
Do(y(x_, 0) = Tanh(x/Sqrt(2));
list(j) = {t, NIntegrate( (1/2)(t-x) Sign(t-x)(y(x,j-1)^2 - 1) y(x, j-1), {x, -10, 10},
AccuracyGoal -> 13, WorkingPrecision -> 13)} // Transpose;
sol(j) = Interpolation(list(j), Method -> "Spline");
y(x_, j) = sol(j)(x)
, {j, m})
```

Basically, I am trivially starting from the exact solution for j=0, and check the convergence at higher order in the iteration (I am doing this to verify whether the code is correct). But from j=3 the numerical solution starts departing from the exact one; so I would like to find a way to correct my code and make the iteration convergent for any higher j, as it should be.

I hope I have clearly stated my problem. Many thanks in advance for your help.