Discretized matrix from the integral kernel function
Recently, I read a paper [1] and then I want to handle the two-dimensional linear integro-differential equation \begin{equation} -\triangle u + q\Big(\frac{\partial u}{\partial x} + \frac{\partial u}{\partial y} \Big) + \gamma\int^{y}_0\int^{x}_0 \mathcal{K}(x, y, s,t) u(s,t)ds dt = g(x, y),\quad\ (x, y) \in \Omega. \end{equation} with the homogeneous Dirichlet boundary condition, where $\Omega = (0, 1)\times (0, 1)$ is the unit square, $\triangle$ is the two-dimensional Laplace operator, $\mathcal{K}(x, y, s, t)$ is the integral kernel function, $g(x, y)$ is a given function, and $q$ and $ \gamma$ are prescribed positive parameters.
After approximating the derivatives by the nine-point centered differences and the integral by the trapezoidal quadrature formula, with the stepsizes being both equal to $h = \frac{1} {\sqrt{m} + 4}$, and normalizing the discretized coefficient matrix and right-hand side by multiplying $h^2$, we obtain a linear system of the form \begin{equation} (T + \zeta B + \nu K)x = b, \end{equation} where $\zeta = qh$ and $\nu = \gamma h^4$ are positive constants. \begin{equation} T = T_1 \otimes I + I \otimes T_1 \in \mathbb{R}^{m\times m}\quad\ \mathrm{and}\quad\ B = B_1 \otimes I + I \otimes B_1 \in \mathbb{R}^{m\times m} \end{equation} are the two-dimensional discrete matrices corresponding to the difference and the integral operators, \begin{equation} T1 = \frac{5}{2} I - \frac{4}{3}(E_1 + E^{T}_1) + \frac{1}{12}(E_2 + E^{T}_2)\quad\ \mathrm{and}\quad\ B_1 = \frac{2}{3}(E_1 - E^{T}_1) - \frac{1}{12}(E_2 - E^{T}_2) \end{equation} are the one-dimensional discrete matrices corresponding to the difference and the integral operators, with $E_1,~E_2\in \mathbb{R}^{\sqrt{m}\times\sqrt{m}}$ being given by, \begin{equation} E_1 = \begin{pmatrix} 0 &1 \ &0 & 1 \ & & 0 & 1 \ & & &\ddots &\ddots \ & & & &0 & 1\ & & & & & 0 \end{pmatrix}\quad\ {\rm and}\quad\ E_2 = \begin{pmatrix} 0 &0 & 1 \ &0 & 0 & 1 \ & &\ddots &\ddots &\ddots \ & & &0 &0 & 1 \ & & & &0 & 0 \ & & & & & 0 \end{pmatrix} \end{equation} and $K = (K_{ij})\in\mathcal{R}^{m\times m}$ is defined by \begin{equation} k_{ij} = x_{i_1}y_{i_2}\mathcal{K}(x_{i_1}, y_{i_2}, x_{j_1}, y_{j_2}), \end{equation} with \begin{equation} i = \sqrt{m}(i_2 - 1) + i_1 \quad\ {\rm and}\quad\ j = \sqrt{m}(j_2 - 1) + j_1,\quad i_1,~i_2,~j_1,~j_2 = 1,2,\ldots,\sqrt{m}. 916