How to solve this coupled integral equation, Mathematica returns the expression without evaluation

Question

ClearAll["Global`*"];
l = 1;
ψ[n_, x_] := 
  Exp[-x^2/(2 l^2)] HermiteH[n, x/l]/Sqrt[2^n Factorial[n] Sqrt[Pi] l];
ee[ξ_, x_, Δ0_, Δ1_] := 
  Sqrt[ξ^2 + 
    0.25 l (ψ[0, x] Δ0 + ψ[2, 
          x] Δ1)^2];
kk[j1_?NumericQ, 
  j2_?NumericQ, ξ_?NumericQ, Δ0_?
   NumericQ, Δ1_?NumericQ] := 
 0.125 NIntegrate[ψ[2 j1, x] ψ[2 j2, x]/
     ee[ξ, x, Δ0, Δ1], {x, -20, 20}]
ξ = 0.1;
NSolve[{kk[0, 
      0, ξ, Δ0, Δ1] Δ0 +
     kk[0, 1, ξ, Δ0, Δ1] \
Δ1 == Δ0, 
  kk[1, 0, ξ, Δ0, Δ1] \
Δ0 + 
    kk[1, 1, ξ, Δ0, Δ1] \
Δ1 == Δ1}, {Δ0, \
Δ1}]

I have checked the above code, I found no problem. It seems Mathematica has not tried to solve it, just return the same expression.

I have tried to solve a simpler equation with only one variable (also have numerical integral), NSolve can handle.

Answer 1

ContourPlot shows the solution:

eqn={kk[0,0, \[Xi], \[CapitalDelta]0, \[CapitalDelta]1] \[CapitalDelta]0 +kk[0, 1, \[Xi], \[CapitalDelta]0, \[CapitalDelta]1] \\[CapitalDelta]1 == \[CapitalDelta]0, kk[1, 0, \[Xi], \[CapitalDelta]0, \[CapitalDelta]1] \\[CapitalDelta]0 + kk[1, 1, \[Xi], \[CapitalDelta]0, \[CapitalDelta]1] \\[CapitalDelta]1 == \[CapitalDelta]1}
pic=ContourPlot[, {[CapitalDelta]0, -.5, .5}, {\[CapitalDelta]1, -.5, .5}];

solutions:

pi=Graphics`Mesh`FindIntersections[pic  ];
Show[{pic,Graphics[{Red,Point[pi]}]}]

altenative

The intersection points might be found with NMinimize instead of NSolve:

J = # . # &[eqn /. Equal -> Subtract]

For example the two positive solutions:

NMinimize[{J, \[CapitalDelta]0 > \[CapitalDelta]1 >0}, {\[CapitalDelta]0, \[CapitalDelta]1}]   
(*{3.13997*10^-21, {\[CapitalDelta]0 -> 0.273414, \[CapitalDelta]1 -> 0.157877}}*)
NMinimize[{J, \[CapitalDelta]1 > \[CapitalDelta]0 >0}, {\[CapitalDelta]0, \[CapitalDelta]1}]
(*{5.69988*10^-18, {\[CapitalDelta]0 -> 0.126083, \[CapitalDelta]1 -> 0.250507}}*)    /