###### Abstract
The eigenfunctions of the following even-order operator  
\[\mathcal{L}_{0}=(-1)^{q}\frac{d^{2q}}{dx^{2q}},\ q\geq 1, \tag{2}\]  
have been investigated by numerous authors (see [58; 70] with references therein) with some Neumann boundary conditions that provide the fastest possible decay of the corresponding expansion coefficients:  
\[\mathcal{L}_{0}\phi(x)=\alpha^{2q}\phi(x),\ x\in[-1,1], \tag{3}\]  
\[\phi^{(r)}(\pm 1)=0,\ r=q,q+1,\ldots,2q-1. \tag{4}\]  
Those eigenfunctions are known as polyharmonic-Neumann eigenfunctions. It is known ([71]) that  
\[\mathcal{L}_{0}\phi(x)=0 \tag{5}\]  
for polynomials of a degree up to \(q-1\). The corresponding orthonormal eigenfunctions are \(\phi_{0,n}\), \(n=0,\ldots,q-1\), where \(\phi_{0,n}\) is the \(n\)-th Legendre polynomial. The remaining eigenvalues are positive and simple. The corresponding system of eigenfunctions  
\[\mathcal{H}_{q}=\{\phi_{0,n}(x):n=0,\ldots,q-1\}\cup\{\phi_{n}(x):n\in\mathbb{ N}\},\ q\geq 1 \tag{6}\]  
forms a dense orthonormal subset of \(L_{2}(-1,1)\). Following [58], we denote the truncated expansion of a function \(f\in L_{2}(-1,1)\) in polyharmonic-Neumann eigenfunctions as  
\[\mathcal{F}_{N}^{q}[f](x):=\sum_{n=0}^{q-1}\hat{f}_{0,n}\phi_{0,n}(x)+\sum_{n =1}^{N}\hat{f}_{n}\phi_{n}(x),\ x\in[-1,1], \tag{7}\]  
where  
\[\hat{f}_{0,n}:=\int_{-1}^{1}f(x)\phi_{0,n}(x)dx,\ \hat{f}_{n}:=\int_{-1}^{1}f(x )\phi_{n}(x)dx. \tag{8}\]  
Approximations by \(\mathcal{F}_{N}^{q}[f](x)\) are better than those based on the classical truncated Fourier series. Let \(H^{r}\) be the Sobolev space with index \(r\). Then, the following estimates are true ([58]):  
1. The coefficients \(\hat{f}_{n}\) of a function \(f\in H^{q+1}(-1,1)\) are \(O(n^{-q-1})\) as \(n\to\infty\) ([58]). 2. If \(f\in H^{q+2}(-1,1)\), the rate of convergence of \(\mathcal{F}_{N}^{q}[f]\) to \(f\) is \(O(N^{-q-1})\) as \(N\to\infty\) in compact subsets of \((-1,1)\) and \(O(N^{-q})\) as \(N\to\infty\) if \(x=\pm 1\). 3. All those results are valid without any periodic conditions, in contrast to the classical truncated Fourier series. The particular case \(q=1\) is of special interest. The corresponding system of eigenfunctions is also known as the modified Fourier basis:  
\[\mathcal{H}_{1}=\{\cos\pi nx:n\in\mathbb{Z}_{+}\}\cup\{\sin\pi(n-\tfrac{1}{2} )x:n\in\mathbb{N}\},\ x\in[-1,1]. \tag{9}\]  
This was proposed by [PERSON] [72] and then explored in papers [73; 74; 75; 76; 77; 78; 70; 73; 74; 76]. Paper [77] considers interpolations on a uniform grid that are exact for the eigenfunctions from the system \(\mathcal{H}_{1}\). It explores the pointwise convergence of the interpolations and proves the convergence rate \(O(N^{-3})\) for compact subsets of \((-1,1)\), which is better than the convergence rate of the corresponding Fourier expansions. Moreover, the paper proves the convergence rate on \([-1,1]\) as \(O(N^{-1})\), the same as for the expansions. The convergence rate of the interpolation in the \(L_{2}(-1,1)\) norm is \(O(N^{-3/2})\), which is again the same as for the corresponding expansions. Hence, the main benefit of applying interpolations is faster convergence at points away from the singularities at \(x=\pm 1\). One of the goals of the current paper is to improve the results of [77] by applying the idea of QP interpolations. We demonstrate how this idea works for interpolations based on the modified Fourier system. We prove some convergence theorems and provide results from numerical experiments. Comparisons with the Fourier expansions and interpolations from [77] confirm our expectations. QP interpolations exhibit improved convergence rates at points away from \(x=\pm 1\). Interestingly, better convergence is also observed over the entire interval, although not necessarily in terms of convergence rate. The next goal is to apply the idea to all polyharmonic-Neumann systems with \(q>1\). While we do not prove convergence theorems for these cases, we provide thorough numerical experiments.