I have an integral that I would like to evaluate over a continuous range of lower limits and plot the result. I can't figure out how to tell Mathematica to do this.

$$Q_c(u) = \frac{2\pi}{c^2} \left(\frac{k_B T_c}{h}\right)^3 \int_a^\infty \frac{u^2}{e^u-1}du$$

Here's my Mathematica code so far.

Q[x_] := x^2/(Exp[x] - 1);
Kc = 2*\[Pi]/c^2*(kb*300/h)^3;
a[vg_] := vg/(kb*300);
z = 1.0;
Integrate[Kc*Q[x], {x, a[z], \[Infinity]}]


It's simple enough to evaluate it for any specific value of $a$ but when I try to make a new function containing the integration Mathematica won't plot the results.

qloop[z_, Q[x_]] := Integrate[Kc*Q[x], {x, a[z], \[Infinity]}];
Plot[qloop[z, Q[x]], {z, .5, 2.5}]

2 Answers 2


You can just plot your integral (whatever Q is, its unit is 1/(m^2 s)) like

c = 3*10^8; kb = 1.38/10^23; T = 300; h = 6.626/10^34; 
NIntegrate[u^2/(Exp[u] - 1), {u, a, Infinity}], {a, 0, 10}, 
AxesLabel -> {"a", "Q"}]

enter image description here


Instead of passing Q[x], pass just z to qloop. Inside qloop, x is set by the integrate itself and passed to Q[x] then. Also change the function to use ?NumericQ.

You did not have values for c and h (I assumes these are speed of light and Planck constant).

Q[x_?NumericQ] := x^2/(Exp[x] - 1);
kb = 1;
c = QuantityMagnitude@UnitConvert[Quantity["SpeedOfLight"]];
h = QuantityMagnitude@UnitConvert[Quantity["PlanckConstant"]];
Kc = 2*π/c^2*(kb*300/h)^3;
a[vg_] := vg/(kb*300);

qloop[z_?NumericQ] := Integrate[Kc*Q[x], {x, a[z], ∞}];
Plot[qloop[z], {z, .5, 2.5}]

Mathematica graphics


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.