Does DCT diagonalize the FD discretisation of the Laplacian with Neumann boundary conditions?

If one has the Poisson problem (assume $$\int_{\Omega} f = 0$$ and $$\int_{\Omega} u = 0$$):

\begin{alignat}{3} \Delta u(x) &= f(x), &\quad&x\in\Omega \\ \partial_nu(x) &= 0, &\quad&x\in\partial\Omega, \end{alignat}

then a standard approach is to discretise it using finite differences. For example in 1D we have:

$$$$(\Delta u)_i \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h^2}.$$$$

Enforcing the Neumann boundary conditions one ends up with a system of the form:

$$$$\begin{bmatrix} -1 & 1 & 0 & 0 & 0 & \ldots & 0\\ 1 & -2 & 1 & 0 & 0 & \ldots & 0\\ 0 & 1 & -2 & 1 & 0 & \ldots & 0\\ &&&\ldots&& \\ 0 & \ldots & 0 & 1 & -2 & 1 & 0\\ 0 & \ldots & 0 & 0 & 1 & -2 & 1\\ 0 & \ldots & 0 & 0 & 0 & 1 & -1\\ \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ \ldots \\ u_n \end{bmatrix} = \begin{bmatrix} f_1 \\ f_2 \\ \ldots \\ f_n \end{bmatrix}$$$$

I have seen (e.g. https://elonen.iki.fi/code/misc-notes/neumann-cosine/) the DCT being suggested as a solution for this problem. If I term the above matrix $$L$$, then as I understand it, the DCT matrix $$V^T$$ is supposed to diagonalize $$L = V\Lambda V^T:$$ \begin{align} Lu&=f \\ V^TLu&=V^Tf \\ V^TV\Lambda V^T u &= V^Tf \\ V^T u &= \Lambda^{-1}V^T f \\ u &= V \Lambda^{-1} V^T f. \end{align} I know for instance that the DFT diagonalizes circulant matrices, however I wasn't aware that the DCT diagonalizes matrices of the above form. I found the following work on the structure of matrices that the discrete cosine and sine transforms diagonalize: "Matrices diagonalized by the discrete cosine and discrete sine transforms". I decided to verify that this is indeed the case. I took as an example: $$$$L = \begin{bmatrix} -1 & 1 & 0 \\ 1 & -2 & 1 \\ 0 & 1 & -1 \end{bmatrix}$$$$ and used the form described in the paper to try and match it: $$$$\begin{bmatrix} t_0 & \sqrt{2}t_1 & t_2 \\ \sqrt{2}t_1 & t_0+t_2 & \sqrt{2}t_1 \\ t_2 & \sqrt{2}t_1 & t_0 \end{bmatrix} \, \textrm{or} \, \begin{bmatrix} t_0 & \sqrt{2}t_1 & \sqrt{2}t_2 \\ \sqrt{2}t_1 & t_0+t_2 & t_1+t_2 \\ \sqrt{2}t_2 & t_1+t_2 & t_0+t_1 \end{bmatrix}.$$$$

However neither of the described forms (corresponding to different types of DCT) can produce the above matrix, as $$t_0=-1, t_2 = 0 \implies t_0+t_2 = -1 \ne -2$$. Am I misunderstanding something?

In fact, since the DCT acts as if the signal had been mirrored, I would expect that it should be able to diagonalize any circulant matrix $$C\in\mathbb{R}^{(n+2m)\times(n+2m)}, \, m\leq n$$ that has been restricted as follows:

$$$$\begin{bmatrix} \boldsymbol{0}_{n\times m} & \boldsymbol{I}_{n} & \boldsymbol{0}_{n\times m} \end{bmatrix} C\begin{bmatrix} \boldsymbol{A}_{m\times n} \\ \boldsymbol{I}_{n} \\ \boldsymbol{B}_{m\times n} \end{bmatrix},$$$$

where $$\boldsymbol{I}_n$$ is the $$n\times n$$ identity matrix and $$\boldsymbol{A}_{m \times n}$$ is the exchange matrix (ones on the antidiagonal) with its top $$(n-m)$$ rows removed, while $$\boldsymbol{B}_{m \times n}$$ is the exchange matrix with it bottom $$(n-m)$$ rows removed. Essentially the above acts as if the vector to be multiplied on the right was reflected on the boundaries for $$m$$ elements, e.g. $$v = (v_1, \ldots, v_n)$$ results in: $$$$u = (v_m, \ldots, v_1, v_1, \ldots, v_n, v_n, \ldots, v_{n-m+1}),$$$$ then the full $$C$$ is applied, and finally the extra elements are thrown away.

The "issue" seems to have been that there is a discrepancy between the used transforms. The DCT transforms discussed in the linked paper are orthogonal, i.e. $$U^{-1}=U^T$$, where the even one corresponds to DCT-I. The orthogonal DCT-I cannot diagonalize the Laplacian (with midpoint reflection, for more details see section 3 in Strang's paper on DCT). However the non-orthogonal version $$V = D U D^{-1}, \, D=\operatorname{diag}(d)$$ for a suitable $$d$$ would diagonalize the Laplacian with midpoint reflection.
Since the matrices considered in the paper have the form $$U \Lambda U^T$$ if I now require that $$L = V\Lambda V^{-1}$$ then this is equivalent to $$L = DUD^{-1}\Lambda D^{-1}U^TD$$, but $$\Sigma = D^{-1}\Lambda D^{-1}$$ is a diagonal matrix also, so the whole thing can be written as $$D(U\Sigma U^T)D$$. Since the structure of $$U\Sigma U^T$$ is given explicitly, all I need to do is multiply by $$D$$: $$$$\begin{bmatrix} (d_1)^2 t_0 & (d_1d_2)\sqrt{2} t_1 & (d_1d_3) t_2 \\ (d_2d_1) \sqrt{2}t_1 & (d_2)^2 (t_0+t_2) & (d_1d_3)\sqrt{2} t_1 \\ (d_3d_1) t_2 & (d_3 d_2) \sqrt{2}t_1 & (d_3)^2 t_0 \end{bmatrix}$$$$
Then for $$d_1 = \frac{1}{\sqrt{2}}, \, d_2 = 1,\, d_3 = \frac{1}{\sqrt{2}}$$ and $$t_0 = -2,\, t_1 = 1,\, t_2=0,$$ I get the desired result. For arbitrary $$n$$ this extends with $$d_1 = 2^{-0.5},\,d_2 = 1,\, \ldots,\, d_{n-1} = 1,\, d_n = 2^{-0.5}$$ and $$t_0 = -2,\, t_1 = 1,\, t_i = 0$$. This rescaling corresponds to the standard non-orthogonal DCT-I transform (if it is additionally multiplied by the scalar $$\sqrt{\frac{2}{n-1}}$$).
Consequently the diagonalization has the form for DCT-I: \begin{align} Lu &= f \\ V^{-1}V \Lambda V^{-1} u &= V^{-1} f \\ u &= V\Lambda^{-1}V^{-1} f. \end{align}
It turns out that the unitary version of DCT-II: $$$$U_{ij} = \frac{C_i}{\sqrt{N}}\cos\left(\frac{\pi i (j+\frac{1}{2})}{N}\right), \, C_1 = 1, \, C_i = \sqrt{2}, \, i\ne 1,$$$$ also diagonalizes the matrix $$L$$ I was considering. The mismatch occurred because the authors of the "Matrices diagonalized by the discrete cosine and discrete sine transforms" paper use a unitary DCT-I matrix in their $$n$$ even case (the notation in that paper is a bit non-standard). I found a very easy to follow write up by Strang on this however. It also shows how other DCT types as well as DST types model different boundary conditions (and mixtures of those).