How to handle this constant-interval integral in ODE?

Question

The problem is to solve a system of nonlinear equations with a definite integral. $Q_i(x), i=1,2,3, x \in [-\text{max},\text{max}]$ \begin{align*} a_iQ_i''(x)-b_i(\vec{Q})-d_i\int_{-\text{max}}^{\text{max}} {\hat K Q_i dx'} = 0 \end{align*}

How can I solve the cases like $d_i=0.01,0.1$ or any other values. It appears to me Laplace transform doesn't help. I failed to turn it to a system of ODEs by differentiation, either. Or I missed something?
So the question is about dealing with the integral. Any suggestions or comments? Thanks in advance.

I use the following code to easily solve the case when all $d_i=0$. The parameters are set to their typical values for my need ($\alpha_1-\alpha_{2,3}>2w_1>0=w_{2,3}$), but one can tune. $b_i=\frac{2}{\alpha_i}(Q_i-w_i)-\frac{Q_i}{|\vec{Q}|}$ and $w_i$ is constant. The form of boundary condition is not to be changed.

The integral kernal $\hat K Q_i=\frac{Q_i(x)-Q_i(x')}{(x-x')^2}$ is singular, although somehow suppressed by the nominator. If not possible to solve this, adding some small cutoff $\frac{Q_i(x)-Q_i(x')}{(x-x')^2+\delta^2}$ is OK.
When all $d_i=0$, it seems that one of $Q_2,Q_3$ stays at zero. I wonder if this holds even when $d_i\neq0$.

ClearAll["Global`*"];
max = 5.0; accur = 15;
w1 = 0.01; w2 = 0; w3 = 0; w = {w1, w2, w3};
α1 = 1.5; α2 = 1.0; α3 = 0.8; α = {\
α1, α2, α3};
Qstart = {w1 - α1/2, 0.0, 0.0}; Qend = {w1 + α1/2, 0.0, 
  0.0};
M1 = 1/α1; M2 = 1/α2; M3 = 1/α3; M = {M1, M2, M3};
dV = {2/α1 (Q1[u] - w1) - Q1[u]/Sqrt[
    Q1[u]^2 + Q2[u]^2 + Q3[u]^2], 
   2/α2 (Q2[u] - w2) - Q2[u]/Sqrt[Q1[u]^2 + Q2[u]^2 + Q3[u]^2],
    2/α3 (Q3[u] - w3) - Q3[u]/Sqrt[
    Q1[u]^2 + Q2[u]^2 + Q3[u]^2]};
s = NDSolve[{M1 Q1''[u] == dV[[1]], M2 Q2''[u] == dV[[2]], 
    M3 Q3''[u] == dV[[3]], Q1[-max] == Qstart[[1]], 
    Q2[-max] == Qstart[[2]], Q3[-max] == Qstart[[3]], 
    Q1[max] == Qend[[1]], Q2[max] == Qend[[2]], 
    Q3[max] == Qend[[3]]}, {Q1, Q2, Q3}, {u, -max, max}, 
   AccuracyGoal -> accur];
QQ = First[{Q1, Q2, Q3} /. s];
Plot[{QQ[[1]][u], QQ[[2]][u], QQ[[3]][u]}, {u, -max, max}, 
 PlotRange -> All, PlotLegends -> Automatic]

Answer 1

The only solution I can think out is to make a step backward i.e. turn to finite difference method. I'll use pdetoae for the generation of difference equation.

V[i_] = (2 (Q[i][x] - w[i]))/α[i] - Q[i][x]/Sqrt@Sum[Q[i][x]^2, {i, 3}];
transrule = lst_List[i_] :> lst[[i]];
pararule =
  {w -> {0.01, 0, 0}, α -> {1.5, 1, 0.8}, d -> {0.01, 0.02, 0.03},
   m -> 1/α, ϵ -> 10^-6};
{Qstart, Qend} = {{w[1] - α[1]/2, 0, 0}, {w[1] + α[1]/2, 0, 0}} //. 
    pararule /. transrule;

Clear@kernel
kernel[Q[i_]@x_] = (Q[i][x] - Q[i][ξ])/((x - ξ)^2 + ϵ) /. pararule;

SetAttributes[int, Listable];
eq = Table[m[i] Q[i]''[x] - V[i] - d[i] int[Q[i]@x] == 0, {i, 3}] //. pararule /. 
   transrule;
max = 5;
points = 100;
difforder = 4;
domain = {-max, max};

{nodes, weights} = Most[NIntegrate`GaussRuleData[points, MachinePrecision]];

midgrid = Rescale[nodes, {0, 1}, domain];

intrule = int[
    qix_] :> -Subtract @@ domain weights.Map[Function[ξ, #] &@kernel[qix], midgrid];

grid = Flatten[{domain // First, midgrid, domain // Last}];

ptoafunc = pdetoae[{Q[1], Q[2], Q[3]}[x], grid, difforder];
del = #[[2 ;; -2]] &;
ae = del /@ ptoafunc[eq] /. intrule;
aebc = Table[
   {Q[i][domain // First] == Qstart[[i]],
    Q[i][domain // Last] == Qend[[i]]}, {i, 3}];

Notice that pdetoae can't handle integral so some extra coding has been made. The code for discretization of integral is modified from this answer.

Last step is to solve the nonlinear equation system with FindRoot:

sollst = With[{initial = -2}, 
   Partition[FindRoot[{ae, aebc}, 
      Flatten[Table[{Q[i][x], initial}, {i, 3}, {x, grid}], 1], 
      MaxIterations -> 200][[All, -1]], points + 2]];

ListLinePlot[{grid, #}\[Transpose] & /@ sollst, PlotRange -> All, Axes -> {False, True}]

Notice there seem to be multiple solutions for the system. With initial = -2:

With initial = 2: