Approximaing Gamma function For $c>1$ and $0<\theta<1$, we wish to approximate (upper bound) following Gamma function:
$$\int_c^{c\theta}x^{-3}e^{-x}dx $$
 A: Since $e^{-x}$ decreases so quickly (and with $\dfrac 1{x^3}$ decreasing quickly for small $x>0$) we may deduce that following integral is approximate (and bounded!) by :
\begin{align}
2\int_{c\theta}^c x^{-3}\,e^{-x}\;dx&\approx \int_{c\,\theta}^\infty 2\,x^{-3}e^{-x}dx,\quad(*)\\
&\approx \left(\frac 1{(c\,\theta)^2}-\frac 1{c\,\theta}\right)e^{-c\,\theta}-\operatorname{Ei}(-c\,\theta)\\
\end{align}
I considered the double of your integral with the sign changed (since $c\theta<c$) making this an upper-bound only for the absolute value of your integral.
$(*)$ integrating by parts twice starting with $\;\displaystyle\frac d{dx}x^{-2}=-\frac 2{x^3}\;$ and from the definition of the exponential integral $\;\displaystyle \operatorname{Ei}(-x)=-\int_{x}^\infty\frac {e^{-t}}t\;dt$
The exponential integral (a common variant is $\operatorname{E1}$ with $\;\operatorname{E1}(x)=-\operatorname{Ei}(-x)\;$ for $x\in\mathbb{R}$) may be easier to evaluate (see A & S and the neat C.F. from DLMF for example).
$$-$$
An exact evaluation is obtained using the difference of two incomplete gamma functions $\displaystyle\Gamma(a,x):=\int_x^\infty e^{-t}\,t^{a-1}\,dt\;$ :
$$\int_{c\,\theta}^c x^{-3}\,e^{-x}\;dx=\Gamma(-2,c\,\theta)-\Gamma(-2,c)$$
Following continued fraction $(6.5.31)$ may help much to evaluate $\Gamma(-2,x)$ (for $x>0$) :
$$\Gamma(-2,x)=\frac{e^{-x}}{x^2}\;\frac 1{x+\cfrac 3{1+\cfrac 1{x+\cfrac 4{1+\cfrac 2{x+\cfrac 5{1+\cfrac 3{x+\cdots}}}}}}}$$
Continued fractions are very interesting for upper (lower) bounds since every additional iteration (fraction) will alternate between too large and too small. We may thus obtain an upper bound for $\Gamma(-2,c\,\theta)$ that we may combine with a lower bound $\Gamma(-2,c)$ (stopping the iteration one step earlier) to obtain an upper bound for the difference. The more terms for the first evaluation the higher the precision!
Hoping this helped,
pari/gp code for incomplete gamma using continued fractions :
ig(a,x,n)=local(r=x+n+1-a);forstep(i=n,1,-1,r=x+(i-a)/(1+i/r));exp(-x)*x^a/r

