The purpose of my code is to sole the Friedmann equation by numerically integrating for a given pair of Omega_M,0 and Omega_Lambda,0. This is my code below and I need help figuring out how to fix my error in my code. This is what I get when I try to run my code.

Traceback (most recent call last):
  File "/home/shakeila/project.py", line 30, in <module>
    a_values, H_values = solve_friedmann(Omega_M0, Omega_Lambda0)
  File "/home/shakeila/project.py", line 18, in solve_friedmann
    result, _ = quad(integrand, 0, a, args=(Omega_m0, Omega_lambda0))
  File "/usr/lib/python3/dist-packages/scipy/integrate/_quadpack_py.py", line 351, in quad
    retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
  File "/usr/lib/python3/dist-packages/scipy/integrate/_quadpack_py.py", line 463, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
TypeError: integrand() takes 2 positional arguments but 3 were given

I try taking out the args= part in my code but I get the same type of error.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

# Constants
H0 = 70.0  # Hubble constant in km/s/Mpc
Omega_r0 = 0.0  # Radiation density parameter (ignored)

# Function to integrate
def integrand(Omega_m0, Omega_lambda0):
    return 1.0 / np.sqrt((Omega_m0 / a) + (Omega_lambda0 * a**2) + 1),

# Function to solve Friedmann Equation
def solve_friedmann(Omega_m0, Omega_lambda0):
    a_values = np.linspace(0.01, 10, 1000)
    a_output = []
    for a in a_values:
        result, _ = quad(integrand, 0, a, args=(Omega_m0, Omega_lambda0))
        a_output.append(result)
    return a_values, np.sqrt(a_output) * H0

# Define Omega_M0 and Omega_Lambda0 pairs
Omega_M0_values = [0.3, 0.5, 0.7]
Omega_Lambda0_values = [0.7, 0.5, 0.3]

# Plotting
plt.figure(figsize=(10, 6))
for Omega_M0 in Omega_M0_values:
    for Omega_Lambda0 in Omega_Lambda0_values:
     a_values, H_values = solve_friedmann(Omega_M0, Omega_Lambda0)
        plt.plot(a_values, H_values, label=f"Ω_M0={Omega_M0}, Ω_Λ0={Omega_Lambda0}")
plt.xlabel('H_0(t-t_0)')
plt.ylabel('a')
plt.title('Scale Factor Evolution for Different Ω_M0 and Ω_Λ0')
plt.legend()
plt.grid(True)                                                     
plt.show()

I tried taking out the args part of the code that did not work. I’m suppose to get a graph of the scale factor for different Omega_M,0 and Omega_Lambda,0.