Performing integral over infinity with numbers in excess of 10^1000

I am attempting to perform on an integral from -infinity to infinity, with an equation containing a term that is raised to the power M. At very high (5000+) values of M, this value of this term can exceed 10^100.

I have this python (3.10) program, which works for M = 1000, but nothing much beyond that.

import numpy as np
from scipy import integrate, special

def integrand(u, M, Z, r):
    sqrt_r = np.sqrt(r)
    sqrt_term = np.sqrt(2 * (1 - r))
    arg1 = (Z + sqrt_r * u) / sqrt_term
    arg2 = (-Z + sqrt_r * u) / sqrt_term
    erf_term = special.erf(arg1) - special.erf(arg2)
    return np.exp(-u**2 / 2) / np.sqrt(2 * np.pi) * erf_term**M

def evaluate_integral(M, Z, r):
    result, _ = integrate.quad(integrand, -np.inf, np.inf, args=(M, Z, r))
    return 1 / (2**M) * result

# Example usage:
M = 1000
Z = 5
r = 0.6



integral_value = evaluate_integral(M, Z, r)

The error I get is:

IntegrationWarning: The occurrence of roundoff error is detected, which prevents 
  the requested tolerance from being achieved.  The error may be 
  underestimated.
  result, _ = integrate.quad(integrand, -np.inf, np.inf, args=(M, Z, r))

and the result for integral_value is nan.

Is there a way I can perform this integral using Python, which would be convenient as it allows me to be included with the rest of my program, or will I be forced to use a different tool?

No, your only problem is that you took the factor 1/2**M outside of the integral. If you put it back inside you will then have

(erf_term/2)**M

and that will not blow up because the difference of two error functions cannot exceed 2, so the new bracketed factor is less than 1.

Assuming that the function is defined properly (which would be surprising given the magnitudes in question), don't integrate it in linear space. Instead use the log equivalent:

def integrand_log(u, M: float, Z: float, r: float):
    sqrt_r = np.sqrt(r)
    sqrt_term = np.sqrt(2*(1 - r))
    arg1 = (sqrt_r*u + Z)/sqrt_term
    arg2 = (sqrt_r*u - Z)/sqrt_term
    erf_term = special.erf(arg1) - special.erf(arg2)
    return (
        M*np.log(erf_term)
        - 0.5 * (u**2 + np.log(2 * np.pi))
    )

This log-identity

gets close to helping you translate back to a linear space integral but is not enough. I'll elaborate in an edit later.