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}. \end{equation*} Here we choose \begin{equation*} \mathcal{K}(x,y,s,t) = \frac{\mu(\sqrt{\psi} + 4)}{xy} \Big[\sqrt{\psi}(y - t) + (x - s) \Big], \end{equation*} with $\psi$ being a positive constant, so that \begin{equation*} K_{ij} = \frac{\mu(\sqrt{\psi} + 4}{\sqrt{m} + 4}\Big[\sqrt{\frac{\psi}{m}}(i - j)+ \Big(1 - \sqrt{\frac{\psi}{m}}\Big) (i_1 - j_1)\Big],\quad\ i,j = 1,2,\ldots,m. \end{equation*} For this problem, the matrix $T$ is symmetric positive definite, and the matrices $B$ and $K$ are skew-symmetric. Besides, both $T$ and $B$ are sparse, but $K$ is dense.
My problems: (1) How to evaluate and assemble the dense matrix $K$ by MATLAB? (There are $i_1,i_2,j_1,j_2,i,j$ make me confused!) (2) How to know that the dense matrix $K$ is skew-symmetric? (3) It seems that we have $i_1 = i_2,~j_1 = j_2$?
