How to calculate joint distribution by integrating over all possible values of model parameters and observations

I was wondering how to calculate the following joint distribution by assuming that $x_i$'s are continuous observations from normal distirbution $N(\mu, \sigma)$ with mean $\mu$ and variance $\sigma$ and $p(\mu, \sigma)$ is the joint distribution of $\mu$ and $\sigma$: $$p(x_1,x_2,..., x_n, x_{n+1}) = \iint \prod_{j=1}^{n+1} N(\mu,\sigma)|_{x_j}} . p(\mu, \sigma) \,d\mu \, d\sigma$$ I know that this can be accomplished by solving the intergration numerically (e.g. with rectangular integration), however, the thing that puzzles me is that the products of $N(\mu,\sigma)|_{x_j}$ would be always zero since the probability of the point $x_i$ in a continuous distribution $N(\mu,\sigma)$ is zero (or I'm wrong?).

In addition, what is the best method to calculate (or estimate) this distribution in python?

$N(\mu,\sigma)|_{x_j}$ is not the probability of a point $x_j$, this is the density of point $x_j$ given the parameters $\mu, \sigma$. The total integral is not the probability either, it is the likelihood of sorts or a density function. The whole thing is in fact never a zero, because the density $N$ of normal distribution is never zero, it can be very small of course.