Fitting using integral

Hello,

I would like to fit using the function
Cv=a9Nk(x/θ)^3 Integral (dx (x^4 e^x)/((e^x-1)^2 )(With limits 0 to θ/x)

Can anyone help me how to write a User defined function for this expression.

I tried to write the integral part alone in the procedure:

#pragma rtGlobals=3 // Use modern global access method and strict wave access.

Function HCinteg(inX)
Variable inX
return (inX^4*e^inX)/((e^inX-1)^2)
End
printf "%.10f|r" Integrate1D(HCinteg, theta/inX, 0)

But the name of the function was not appearing in the Analysis/curve fitting menu.
I have attached the image of the function clearly.

Thank you very much