Solving double integral numerically in matlab In the paper "The fractional Laplacian operator on bounded domains as a special case of the nonlocal diffusion operator". Where the author has solved a fractional laplacian equation on bounded domain as a non-local diffusion equation.
I am trying to implement the finite element approximation of the one dimensional problem(please refer to page 14 of the above mentioned paper) in matlab.
I am using the following definition of $\phi_k$ as it is mentioned in the paper that $\phi$ is a $hat\;function$
\begin{equation}
\phi_{k}(x)
=
\begin{cases}
{\large{x-x_{k-1} \over x_k\,-x_{k-1}}}
&
\mbox{ if } x \in [x_{k-1},x_k],
\\[3mm]
{\large{x_{k+1}\,-x \over x_{k+1}\,-x_k}}
&
\mbox{ if } x \in [x_k,x_{k+1}],
\\[3mm]
0 & \mbox{ otherwise},
\end{cases}
\end{equation}
$\Omega=(-1,1)$ and $\Omega_I=(-1-\lambda,-1) \cup (1,1+\lambda)$ so that $\Omega\cup\Omega_I=(-1-\lambda,1+\lambda)$
For the integers K,N we define the partition of $\overline{\Omega\cup\Omega_I}=[-1-\lambda,1+\lambda]$ as,
\begin{equation}
-1-\lambda=x_{-K}<...<x_{-1}<0=x_{0}<x_{1}<...<x_{N-1}<x_{N}=1<x_{N+1}<...<x_{N+K}=1+\lambda
\end{equation}
Finally the equations that we have to solve to get the solution $\tilde{u_N}=\sum_{i=-K}^{K+N}U_j\phi_j(x)$ for some coefficients $U_j$ is:
\begin{equation}
\frac{c_{n,s}}{2}\sum_{j=1}^{N-1}U_{j}\intop_{-1-\lambda}^{1+\lambda}\intop_{x-\lambda}^{x+\lambda}\frac{\phi_{j}(y)-\phi_{j}(x)}{|y-x|^{1+2s}}(\phi_{i}(y)-\phi_{i}(x))dydx=\intop_{-1}^{1}f\phi_{i}dx
\end{equation}
Where $i=1,...,N-1$. 
I need pointers in order to simplify and solve the LHS double integral in matlab.
It is written in the paper(page 15) that I should use four point gauss quadrature for inner integral and quadgk.m function for outer integral, but since the limits of the inner integral are in terms of x how can I apply four point gauss quadrature on it?
I also wrote to the author and she said: 
Solving the integrals in the matrix of the system is the difficult part in solving numerically nonlocal problems. The idea is to implement 2 different quadrature rules for the inner and the outer integral. 
Inner integral: this one goes from x-epsilon to x+epsilon. Since we want to avoid the singularity at y=x we split the integral in (x-epsilon,x) and (x,x+epsilon). Then, take e.g. (x-epsilon,x) and split it according to the elements in your partition. You'll have something like
(x-epsilon, x_j)
(x_j, x_{j+1})
(x_{j+1}, x_{j+2})
(x_{j+2}, x)
for x_j, x_{j+1} and x_{j+2} in (x-epsilon, x).
In these intervals you can use a Gauss quadrature rule (unless you want to have a very very accurate solution, then you should use an adaptive quadrature rule). 
Outer integral: here you need a more accurate rule. You have to integrate over each element of your partition, i.e. (x_i, x_{i+1}) and you should use an adaptive quadrature rule (like quadgk in Matlab). The integrand function is the result of the inner integration.
How should I go about implementing it? 
Any help will be appreciated. Thanks.
 A: Denote 
$$A_{i, j} = \intop_{-1-\lambda}^{1+\lambda}\intop_{x-\lambda}^{x+\lambda}\frac{\phi_{j}(y)-\phi_{j}(x)}{|y-x|^{1+2s}}(\phi_{i}(y)-\phi_{i}(x)) \, dy \, dx.$$
Observe that $A_{i, j}$ can be written as the a single integral:
$$A_{i, j} = \int_{-1 - \lambda}^{1 + \lambda} F(x) \, dx,$$
where
$$F(x) = \intop_{x-\lambda}^{x+\lambda}\frac{\phi_{j}(y)-\phi_{j}(x)}{|y-x|^{1+2s}}(\phi_{i}(y)-\phi_{i}(x)) \, dy.$$
Therefore $A_{i, j}$ can be implemented in MATLAB using quadgk e.g. as Aij = quadgk(@myfun, -1 - lambda, 1 + lambda), where myfun.m is a function which returns the value of $F(x)$. quadgk uses an adaptive quadrature rule, so it should be accurate (assuming that your calculation of $F(x)$ in myfun.m is accurate).
As stated by the author, the difficulty is then in the evaluation of $F(x)$ (alternatively, in the myfun.m file). First you must split the integral into two parts -- the reason why is that most numerical quadrature rules, even if they are designed specifically for the type of singularity under consideration, require the singularity to be at the end points of the interval, otherwise the numerical answer you compute might be completely bogus.
Therefore, split the integral into two intervals: $[x - \lambda, x - \epsilon]$ and $[x + \epsilon, x + \lambda]$ (where $\epsilon$ is a small parameter satisfying $0 < \epsilon < \lambda$ that you would assign at the beginning of the calculations). There may be some small error due to neglecting the contribution to the integral on the interval $(x - \epsilon, x + \epsilon)$, but this small error is a necessary evil if there is no known quadrature rule that handles the singularity $|y - x|^{1 + 2s}$.
Next, recall that the FEM basis functions you are using are linear hat functions: they are linear (or zero) on the intervals $[x_i, x_{i + 1}]$, $i = -K, -K + 1, \ldots, N + K - 1$. Therefore you should further split up the intervals $[x - \lambda, x - \epsilon]$ and $[x + \epsilon, x + \lambda]$ according to the grid:
$$[x - \lambda, x_j], [x_j, x_{j + 1}], \ldots, [x_{j + k}, x - \epsilon]$$
and
$$[x + \epsilon, x_l], [x_l, x_{l + 1}], \ldots, [x_{l + m}, x + \lambda].$$
The notation is cumbersome, but what I mean is this: suppose that $x_j < x_{j + 1}, \ldots x_{j + k}$ are all of the grid points contained in the interval $[x - \lambda, \ldots, x - \epsilon]$, and that $x_l < x_{l + 1} < \ldots < x_{l + m}$ are all of the grid points contained in the interval $[x - \epsilon, x + \lambda]$. Splitting up the two integrals on $[x - \lambda, x - \epsilon], [x + \epsilon, x + \lambda]$ onto these intervals allows you to use a Gauss quadrature rule on each subinterval (as mentioned by the author) which should be fairly accurate, because of the nature of $\phi_i, \phi_j$ (either zero or a linear function on each subinterval)
