I have attempted to approximate this Gaussian:
$$ I =\int_{0}^{\lambda}dx\left(r+x\right)\exp\left(-\rho\left(ax^{1/2}+bx^{3/2}\right)\right)\ $$
using
$$ I =\int_{0}^{\lambda}dx\left(r+x\right)e^{-\rho\left(ax^{1/2}\right)}\left(1-b\rho x^{3/2}\right)\ $$
where the exponent
$$ ax^{1/2}+bx^{3/2}\ $$
is only really an approximation to a much more complicated function, valid for small $x$ (hence the approximation in the second equation).
The question is, I may well be loosing too much accuracy by making this approximation. Can anyone think of a clever or simply better way to approximate $I$? Or is my approximation the best I can do, given the conditions?
If you have only a few of these integrals to do, or if you want to evaluate a few "benchmarks" to test a self-written software implementation, then I'd recommend an interactive CAS with numeric capabilities such as Maxima or Sage.
You would only be able to develop a competitive numerical integration routine under special circumstances, e.g. possessing knowledge of the integrand that can be exploited for improved approximation. The theory of numerical integration (quadrature) is not deep and is often treated in introductory numerical analysis courses. But coding and testing of adaptive quadrature routines has reached a high state of the art in many off-the-shelf libraries. Why reinvent the wheel?
What does make sense is to take a close look at the parameters of your integral to see if some of them can be removed or combined. For example, the parameter $r$ appears in a relatively trivial way:
$$ I =\int_{0}^{\lambda} \exp\left(-\rho\left(ax^{1/2}+bx^{3/2}\right)\right) \left(r+x\right) dx $$ $$ = r \int_{0}^{\lambda} \exp\left(-\rho\left(ax^{1/2}+bx^{3/2}\right)\right) dx + \int_{0}^{\lambda} \exp\left(-\rho\left(ax^{1/2}+bx^{3/2}\right)\right) x dx $$
Furthermore $\rho$ can obviously be combined with $a$ and $b$, so that using $A = \rho a$ and $B = \rho b$ the integrals can be rewritten in the forms:
$$ \int_{0}^{\lambda} \exp\left(-Ax^{1/2}-Bx^{3/2}\right) x^k dx $$
for $k=0,1$, and both integrals have three parameters: $\lambda, A, B$. We can go further and make a substitution $u = \lambda x$ to remove the variable upper limit of integration, but since you've indicated this was only a simplified version of the integrals, I won't try to carry the discussion (of parameter reduction) to a final conclusion.
Some open source (GPL licensed) numerical integration routines are available in C, derived from Fortran routines written for
QUADPACK(in the public domain, available in Netlib).In some situations, perhaps where the number of parameters can be reduced to one or two, you might be able to fit a polynomial or rational approximation to a finite set of values produced with high-precision adaptive quadrature routines, and thus replace the overhead of calling the quadrature routines themselves with calls to these efficient approximations.