So, here's the deal. I need to solve the paraxial diffraction integral for resonators by using a Gaussian quadrature in order to obtain the modes. Typically, this is done by assuming radial symmetry and transforming it to a 1D integral. However, I can't assume symmetry.

For 1D (Assuming symmetry):

\[\lambda u_{lp}(r_2)=\int_0^{a_1} K_l(r_1, r_2)u_{lp}(r_1)dr_1\]

Then, we take \(\textbf{r}=[r_1, r_2, ..., r_N]\) and \(\textbf{w}=[w_1, w_2,...,w_N]\) to be the abscissa and the weights for the Gaussian quadrature. If we consider equal domains for the variables \(x_1\) and \(x_2\) and \(\textbf{u}=[u_1, u_2, ..., u_N]\) the column vector of the eigenfunction so that \(u_N = u(x_N)\), we can rewrite the integral as

\[\lambda \textbf{u}=(\textbf{H}*\textbf{W})*\textbf{u}\]

Where **H **is a *N x N* matrix where \(H_{m,n}=H(r_m, r_n)\) and **W** is diagonal matrix with the elements of **w**.

Hence, we only need to find the eigenvalues and eigenvectors for the equation before and it's done.

However, things aren't that easy for 2D. The integral equation reads as:

\[\lambda u(x_2, y_2)=\int_0^{a_1} K(x_1, x_2, y_1, y_2)u(x_1, y_1)dx_1dy_1\]

I'm not sure how to transform that integral into the matrix equation and what does \(\textbf{H}\) should be in this case, since it needs to be a 4D array? \(H_{xm, xn, ym, yn}=H(x_m, x_n, y_m, y_n)\)?