I am trying to do the integration in this way. This is just an example but I want to do something like this, the output of the first integration will be in the second integration. I hope it is clear from this. Is there any way to do it?

f[r_] := Sqrt[r^2 + r^3];
p71[r_]:=r^2
a1 = 5.93999103;
g91 = { };
While[a1 <= 20,
p31 = NIntegrate[f[r], {r, 5.93999103, a1}];
p41 = p71[a1]*p31;

AppendTo[g91, p41];
a1 += 0.1];

p318[nu_] :=
NIntegrate[((4*nu)/((Exp[(nu)/(g91)^(1/4)]) - 1)), {r,
5.939991023, \[Infinity]}];
Plot[p318[nu], {nu, 0.001, 10}]

• p71 is undefined Aug 25, 2022 at 14:21
• @yarchik sorry for that. I edited it.
– AAA
Aug 25, 2022 at 14:24
• When I paste your code into a fresh notebook and change only the last line to Table[p318[nu],{nu,1,10}] I get a flood of error messages, most of which look like they are from dividing by zero. Is that what you see? Can you track down why you are getting those?
– Bill
Aug 25, 2022 at 15:06
• @Bill that is only because the first element in g91 is zero. If we ignore that and do the integration then still I am not getting anything. so I do not know if this is the correct way to do it. Since g91 is an array which has elements at every r and I want to integrate new function.
– AAA
Aug 25, 2022 at 15:09
• g91is a numerical list of numbers and occurs in your last integral Exp[(nu)/(g91)^(1/4)] . What could be the purpose of such a construct? Aug 25, 2022 at 15:11

If I correctly understand your calculations, you are looking for a function a0= 5.93999103; F[a1]=a1^2 Integrate[f[r],{r,a0,a1}]?

Instead of While try NDSolve:

a0= 5.93999103;
f[r_] := Sqrt[r^2 + r^3];
D[a1^2 Integrate[f[r],{r,a0,a1}],a1]; (* F'[a1]*)


F'[a1] depends on F[a1]

ode=Derivative[F][a1] == a1^2 f[a1] + (2 F[a1])/a1
g91 = NDSolveValue[{ode, F[a0] == 0}, F, {a1, a0, 20}] (* function g91[a1]*)
Plot[g91[a1], {a1, a0, 20}, AxesLabel -> {a1, "g91[a1]"}] Knowing g91 it's possible (hopefully) to evaluate the second integral.

But second integral doesn't converge!