# How to numerically solve a differential equation which has integration in it?

I am trying to solve the following differential equation which has integration in it. It is not giving any results. Can anyone help?

F[t_, k_, \[Chi]_, x_] :=
t^k*Sqrt[1 + (1817/2) t]/(1 + Exp[t - (\[Chi]/(x*10^7))]);
ODE = \[Chi]'[x] +
4*\[Pi]*x (((3/(4*\[Pi]))*
HeavisideTheta[
1.81712 - x]) - ((Sqrt[2]/\[Pi]^2)*1817^(3/
2)*(Integrate[
F[t, 1/2, \[Chi][x], x], {t,
0, \[Infinity]}] + ((1817)*(Integrate[
F[t, 3/2, \[Chi][x], x], {t, 0, \[Infinity]}]))))) == 0;
sol = NDSolve[{ODE, \[Chi][0] == 0}, \[Chi], {x, 0, 1.81712}];
Plot[Evaluate[\[Chi][x] /. First@sol], {x, 0, 1.81712}]

• The integral does not seem to be intergrable. !Mathematica graphics Commented Feb 24, 2023 at 12:41
• now i have edited it. It is integrable. how do you generally solve such differential equations?
– AAA
Commented Feb 24, 2023 at 12:49
• Could you show how you integrated it? It is still not integrable. !Mathematica graphics without being able to integrate it, no point of going to the next step and use NDSolve. Commented Feb 24, 2023 at 12:52
• that is the problem actually. I need to integrate it at each step since the differential equation is d\chi/dx.. So at each \chi and x , I need to find integration.
– AAA
Commented Feb 24, 2023 at 12:53
• But you are using the thing you are trying to solve for, which is $\chi(x)$ in the integrand. So it is unknown. It is what you are solving for in the differential equation. So this will not work. it is like a catch 22. to solve for it, you need to evaluate an integral which uses the solution. How could this work? At least I do not know. May be someone else will have an idea. Commented Feb 24, 2023 at 12:57

Some progress can be made as follows. First, rescale X by 10^7 to avoid large values of X (at least at first). Then, because X and x appear only as a ratio in the integrands, define the functions,

f1[y_?NumericQ] := NIntegrate[Sqrt[t + (1817/2) t^2]/(1 + Exp[t - y]),
{t, 0, Infinity]}, MaxRecursion -> 100]
f2[y_?NumericQ] := NIntegrate[t Sqrt[t + (1817/2) t^2]/(1 + Exp[t - y]),
{t, 0, Infinity}, MaxRecursion -> 100]


Plots of these two functions are

and the ODE is solved by

xm = Rationalize[.001438, 0];
sol = NDSolveValue[{X'[x] + 4*π*x 10^-7 (3/(4*π) -
(Sqrt[2]/π^2*1817^(3/2)*(f1[X[x]/x] + (1817*(f2[X[x]/x]))))) == 0,
X[10^-8] == 0}, X[x], {x, 10^-8, xm}, WorkingPrecision -> 30];


The lower bound of integration is set to 10^-8 to avoid division by zero at x = 0. The upper bound is set to 0.001438, because the solution blows up at about that value of x, as we see in to following plot.

Plot[sol, {x, 10^-8, xm}, PlotRange -> All, AxesLabel -> {x, X},
ImageSize -> Large, LabelStyle -> {15, Bold, Black}]


The explosive growth shown follows from the ODE and is not a numerical artifact, although accuracy deteriorates with this rapid growth.

Approximate expressions for f1 and f2

For t large relative to 2/1817, f1 can be approximated by

Integrate[Sqrt[1817/2] t/(1 + Exp[t - y]), {t, 0, Infinity}]
(* -Sqrt[1817/2] PolyLog[2, -E^y] *)


and similarly for f2

Integrate[Sqrt[1817/2] t^2/(1 + Exp[t - y]), {t, 0, Infinity}]
(* -2 Sqrt[1817/2] PolyLog[3, -E^y] *)


Unfortunately, this approximation does not help much in solving the ODE.

• how did you plot X vs x. I tried this "Plot[X[x] /. sol, {x, 10^-8, xm}]" and got this error "ReplaceAll::reps: {6.02368*10^-7} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing."
– AAA
Commented Feb 25, 2023 at 5:30
• @AKU I added the code that produces the final plot to my answer. Note that I used NDSolveValue, not NDSolve, to produce sol. Commented Feb 25, 2023 at 14:56

Not an answer, perhaps a starting point...

You could try to find an iterative solution (Picard iteration) using NestList starting with \[Chi][x]=0

sol = NestList[
NDSolveValue[{\[Chi]'[x] +
4*\[Pi]*x (((3/(4*\[Pi])) ) - ((Sqrt[2]/\[Pi]^2)*1817^(3/
2)*(NIntegrate[
F[t, 1/2, #[x], x], {t,
0, \[Infinity]}] + ((1817)*(NIntegrate[
F[t, 3/2, #[x], x], {t, 0, \[Infinity]}]))))) ==
0, \[Chi][0] == 0}, \[Chi], {x, 0, 1.81712}] &, 0 &, 2]


Unfortunately only the first iteration evaluates!

• I tried but as you said, it is not giving the desired result. BTW what is &,2 used for?
– AAA
Commented Feb 24, 2023 at 14:05
• & defines a so called pure function, that' s a function without argument. Where does your ode come from? Perhaps it's possible to modify the expression to get a convergent solution? Commented Feb 24, 2023 at 14:36
• Oh okay. I have this analytical expression only. I do not know if we can do any modification to function F since this is the direct analytical expression we have.
– AAA
Commented Feb 24, 2023 at 18:09
• Where do the integral-parts in your ode come from? Commented Feb 24, 2023 at 18:30
• this is just an expression I am getting like this and it is present like this only.
– AAA
Commented Feb 24, 2023 at 18:39

Perhaps StreamPlothelps to understand the problem little bit better

F[t_, k_, \[Chi]_, x_] := t^k*Sqrt[1 + (1817/2) t]/(1 + Exp[t -(\[Chi]/(x*10^7))]);
int[z_?NumericQ, k_?NumericQ] := NIntegrate[F[t, k, z], {t, 0, \[Infinity]}]

StreamPlot[{1, -4*\[Pi]*x (((3/(4*\[Pi]))*HeavisideTheta[1.81712 - x]) -
((Sqrt[2]/\[Pi]^2)*1817^(3/2)*( int[ \[Chi] /x,1/2] + ((1817)*( int[ \[Chi] /x, 3/2] )))))}
, {x, 0.001, 1.81712 }, {\[Chi], -1, 1}, FrameLabel -> {x, \[Chi]}]


It looks like there is no "stable" solution! Phaseplot also confirms @bbgodfrey's explosive growth observation!