So I have a series of non-linear equations which I wish to solve as fast as possible, to illustrate for the case of $n = 4$, I have the following equations:
\begin{gather*} c_0\exp(\lambda_0-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = 1 \\ c_1\exp(\lambda_0+\lambda_1-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_1 \\ c_2\exp(\lambda_0+\lambda_2-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_2 \\ c_3\exp(\lambda_0+\lambda_3-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_3 \\ c_4\exp(\lambda_0+\lambda_4-1)+f\exp(\lambda_0 + \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 - 1) = P_4 \end{gather*} where \begin{gather*} f = \int_{d_4}^{\infty} \int_{d_3}^{\infty} \int_{d_2}^{\infty} \int_{d_1}^{\infty} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ c_0 = \int_{-\infty}^{d_4} \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ c_1 = \int_{-\infty}^{d_4} \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{d_1}^{\infty} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ c_2 = \int_{-\infty}^{d_4} \int_{-\infty}^{d_3}\int_{d_2}^{\infty} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ c_3 = \int_{-\infty}^{d_4} \int_{d_3}^{\infty}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ c_4 = \int_{d_4}^{\infty} \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, x_4) dx_1 dx_2 dx_3 dx_4 \\ \end{gather*}
The definitions are as follows:
$P_1, P_2, P_3, P_4$ are known values.
$d_1, d_2, d_3, d_4$ are known values.
$q(x_1, x_2, x_3, x_4)$ is a multivariate normal distribution with mean vector $0$ and a known variance-covariance matrix.
Thus $f, c_0, c_1, c_2, c_3, c_4$ are all simply real numbers.
For the general case I have:
\begin{gather*} c_0\exp(\lambda_0-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = 1 \\ c_1\exp(\lambda_0+\lambda_1-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_1 \\ c_2\exp(\lambda_0+\lambda_2-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_2 \\ c_3\exp(\lambda_0+\lambda_3-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_3 \\ \vdots \\ c_n\exp(\lambda_0+\lambda_n-1)+f\exp\left(\sum_{i=0}^n \lambda_i - 1\right) = P_n \end{gather*} where \begin{gather*} f = \int_{d_n}^{\infty} \cdots \int_{d_3}^{\infty} \int_{d_2}^{\infty} \int_{d_1}^{\infty} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ c_0 = \int_{-\infty}^{d_n} \cdots \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ c_1 = \int_{-\infty}^{d_n} \cdots \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{d_1}^{\infty} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ c_2 = \int_{-\infty}^{d_n} \cdots \int_{-\infty}^{d_3}\int_{d_2}^{\infty} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ c_3 = \int_{-\infty}^{d_n} \cdots \int_{d_3}^{\infty}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ \vdots \\ c_n = \int_{d_n}^{\infty} \cdots \int_{-\infty}^{d_3}\int_{-\infty}^{d_2} \int_{-\infty}^{d_1} q(x_1, x_2, x_3, \cdots, x_n) dx_1 dx_2 dx_3 \cdots dx_n \\ \end{gather*}
The definitions are as follows:
$P_1, P_2, P_3, \cdots, P_n$ are known values.
$d_1, d_2, d_3, \cdots, d_n$ are known values.
$q(x_1, x_2, x_3, \cdots, x_n)$ is a multivariate normal distribution with mean vector $0$ and a known variance-covariance matrix.
Thus $f, c_0, c_1, c_2, c_3, \cdots, c_n$ are all simply real numbers.
At most $n$ will be 35, so there will be 36 equations to solve simultaneously. What is the best method to solve this as quickly as possible using whatever programming language possible? Also I will need to change parameter inputs for $P_i$, $d_i$ and the variance covariance matrix and solve for a new set of $\lambda$'s, I will need to do this at least 2000 times, which means effectively I will need to solve a set of 35 equations for 2000 times.
