Getting plots form differential equation with two discrete parameters

Question

I am solving the following differential equation for different values of a and n:

$\qquad y'=- \lambda (y^2-Y_{eq}^2)$

where

$\qquad Y_{eq}(x)=1/\pi^2 \int_0^\infty y^2/[\exp(\sqrt{x^2+y^2})+1]dy$

and

$\qquad \lambda=10^{9+a} x^{-n-2}$

for $a\in\{-2,0,2\}$ and $n\in\{0,1,2\}$. Initial condition is $y(1)=0.001$.

My code is the following:

Yeq[x_] := Integrate[1/π^2 y^2/(Exp[Sqrt[x^2 + y^2]] + 1), {y, 0, Infinity}];

Table[
  NDSolve[
    {y'[x] == -10^(9 + a) x^(-n - 2) (y[x]^2 - Yeq[x]^2), 
     y[1] == 0.001}, 
    y, {x, 1, 100}], 
  {n, {0, 1, 2}}, {a, {-2, 0, 2}}];

How can I get plotted solutions?

Answer 1

If Yeq[x_] is not associated with y[x], then equation is ODE and solution is

Y[x_?NumericQ] := 
  NIntegrate[
   1/\[Pi]^2 s^2/(Exp[Sqrt[x^2 + s^2]] + 1), {s, 0, Infinity}];


Do[sol[n, a] = 
    NDSolve[{y'[x] == -10^(9 + a) x^(-n - 2) (y[x]^2 - Y[x]^2), 
      y[1] == 0.001}, y, {x, 1, 18}];, {n, {0, 1, 2}}, {a, {-2, 0, 
    2}}];

Visualisation in logarithmic scale. We don't need so large interval since solution decreases to zero very fast for a given set of parameters

Table[LogLogPlot[y[x] /. sol[n, a][[1]], {x, 1, 18}, 
  PlotLabel -> {Row[{"n = ", n}], Row[{"a = ", a}]}, PlotRange -> All,
   Frame -> True, Axes -> False], {n, {0, 1, 2}}, {a, {-2, 0, 2}}]