Solving system of integral equations using fsolve and quad in python scipy

I am trying to solve the following system of integral equations with unknowns theta1, theta2, theta3:

where Phi and phi are the cdf and pdf of the standard normal distribution respectively by using scipy's fsolve and integrate. Here is my code:

import numpy as np
import math
import scipy.integrate as integrate
from scipy import integrate
from scipy.stats import norm
from scipy.optimize import fsolve


def function(xi, thetai, thetaj, thetak):
  return (1 - norm.cdf(xi - thetaj)) * (1 - norm.cdf(xi - thetak)) * norm.pdf(xi - thetai)

def pi_i(thetai, thetaj, thetak):
  return integrate.quad(function, -np.inf, np.inf)[0]

def equations(p):
    t1, t2, t3 = p
    return (pi_i(t1,t2,t3) - 0.5, pi_i(t2,t1,t3) - 0.3, pi_i(t3,t1,t2) - 0.2)

t1, t2, t3 =  fsolve(equations, (1, 1, 1))

print(equations((t1, t2, t3)))

However, when I run my code, the following error pops up:

TypeError: function() missing 3 required positional arguments: 'thetai', 'thetaj', and 'thetak'

That's because you need to include the arguments for function function when calling integrate.quad using the args param

def pi_i(thetai, thetaj, thetak):
  return integrate.quad(function, -np.inf, np.inf, args=(thetai, thetaj, thetak))[0]

Note: You also need to fix the issue with theta variable in the function function. It's not provided in the function parameters and not a global variable so I'm not sure where it's from