The integro-differential equation can be solved with only a modest amount of code, as follows: First note that the integrand can be simplified because, as stated in the question, f[r] == g[r]
.
int = (Sin[(y[t] - y[tp])/2] - Sin[(y[t] + y[tp])/2]) f[t - tp];
Simplify[int]
(* -2 Cos[y[t]/2] f[t - tp] Sin[y[tp]/2] *)
and -2 Cos[y[t]/2]
can be moved from within the integral. Additionally, because y[tp]
vanishes for tp <= 0
, the lower bound on the integral can be increased from -Infinity
to 0
. Consequently, the integro-differential equation reduces to
D[y[t], t] == i0 + 2 Cos[y[t]/2]
NIntegrate[Sin[y[tp]/2] (t - tp)/(1 + (t - tp)^2)^2, {tp, 0, t}]
Unfortunately, NDSolve
is unable to handle such equations. It turns out, though, that the equation can be solved iteratively. For instance, with i0 = 1
,
Clear[s, h];
n = 14; tmax = 50; i0 = 1;
s[0] = FunctionInterpolation[i0 t, {t, 0, tmax}];
Do[
h[t_?NumericQ] := NIntegrate[Sin[s[i - 1][tp]/2]
(t - tp)/(1 + (t - tp)^2)^2, {tp, 0, t}];
s[i] = NDSolveValue[{D[y[t], t] == i0 + 2 Cos[y[t]/2] h[t], y[0] == 0},
y, {t, 0, tmax}], {i, 1, n}];
Plot[Evaluate@Array[s[#][t] &, n + 1, 0], {t, 0, tmax}, ImageSize -> Large,
AxesLabel -> {t, y}, LabelStyle -> {15, Bold, Black}]

Successive iterations lead to progressively lower curves, with good convergence obtained with ten iterations. Values at t = tmax
illustrate the convergence.
Array[s[#][tmax] &, n + 1, 0]
(* 50., 48.5208, 47.4234, 45.481, 43.3451, 42.2148, 41.175, 39.5509,
37.8209, 36.9781, 36.6383, 36.4922, 36.4351, 36.4158, 36.41} *)
The computation takes about an hour.
Addendum: Asymptotically constant solutions.
Suppose that y[t]
is approximately constant for large t
. Then, Sin[y[t]/2]
can be removed from within the integral, leaving
Integrate[(t - tp)/(1 + (t - tp)^2)^2, {tp, 0, t}, Assumptions -> t > 0]
{* t^2/(2 + 2 t^2) *}
which reduces to 1/2
for large t
, and the integro-differential equation collapses to i0 == -Sin[y/2]/2
. Hence, asymptotically flat solutions exist if and only if -1/2 < i0 < 1/2
, and the asymptotic values should be, for i0 = {1/10, 2/10, 3/10, 4/10, 5/10}
,
N@Table[Min@SolveValues[{Sin[x]/2 == -j, Pi < x < 2 Pi}, x], {j, 1/10, 5/10, 1/10}]
(* {3.34295, 3.55311, 3.78509, 4.06889, 4.71239} *)
The earlier code block gives corresponding converged solutions at tmax = 50
of
(* {3.34306, 3.55333, 3.78547, 4.06955, 4.61028} *)
Evidently, the i0 = 1/2
solution is not yet in the asymptotic regime at tmax = 50
. Corresponding converged solutions are
