Boys function for Gaussian integrals in ab-initio calculations
I am not aware of any existing Fortran code for direct numerical quadrature of this problem, but it is worth pointing out that Mathematica can perform this integral symbolically:

Integrate[t^(2 n) Exp[-x t^2], {t, 0, 1}]

(* 1/2 x^(-(1/2) - n) (Gamma[1/2 + n] - Gamma[1/2 + n, x]) *)


where the $\Gamma$ function can be computed by exponentiating easy-to-find $\ln \Gamma$ functions, often called lngamma().  (One of them is the Euler function $\Gamma(z)$ and the other incomplete $\Gamma(a,z)$ function).

As a check:

With[
  {x = 1, n = 2},
  1/2 x^(-(1/2) - n) (Gamma[1/2 + n] - Gamma[1/2 + n, x])
  ] // N

(* 0.100269  *)


and numerical quadrature yields:

With[
 {x = 1, n = 2},
 NIntegrate[t^(2 n) Exp[-x t^2], {t, 0, 1}]
 ]

(* 0.100269  *)


One has to be a little bit careful, as the $\Gamma$ function is closely related to the factorial, and it can explode as the arguments get too large.  I would experiment and see whether you can get by with this analytic solution, which essentially takes you from numerical quadrature into evaluation of special functions.

The $\Gamma$ function approach is faster by over two orders of magnitude  than the quadrature approach.  This is critical when evaluating molecular integrals.  I have seen implementations where the Boys function is actually interpolated rather than evaluated, and there is discussion about the error in interpolation being almost insignificant if suitable interpolation basis functions are used.