I want to find a general formula for the trace of the following $N\times N$ matrix raised to the power of $d$, where $d \in \mathbb{N}$. $$ \begin{bmatrix} N-2 & 1 & 0 & \cdots & 0 & 0 \\ 1 & N-3 & 1 & \cdots & 0 & 0 \\ 0 & 1 & N-3 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & N-3 & 1 \\ 0 & 0 & 0 & \cdots & 1 & N-2 \end{bmatrix}^d $$
For example when $N=4$ and $d=2$, it will look like this $$ \begin{bmatrix} 2 & 1 & 0 & 0 \\ 1 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 2 \\ \end{bmatrix}^2 = \begin{bmatrix} 5 & 3 & 1 & 0 \\ 3 & 3 & 2 & 1 \\ 1 & 2 & 3 & 3 \\ 0 & 1 & 3 & 5 \\ \end{bmatrix}, $$ so the value of trace for $N=4$ and $d=2$ is $5+3+3+5=16$.
So far I found the solutions for $N=2$ and $N=3$ and $N=4$, but I can't find any patterns from it. $$ N=2 : a_d=(-1)^d+1 $$ $$ N=3 : a_d=2^d+(-1)^d+1 $$ $$ N=4 : a_d=3^d+(1+\sqrt2)^d+(1-\sqrt2)^d+1 $$
The matrix $$ A = \pmatrix{ -1 & 1 & 0 & \cdots & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 \\ 0 & 1 & -2 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -2 & 1 \\ 0 & 0 & 0 & \cdots & 1 & -1 } $$ corresponds to the discretized second derivative operator with Neumann boundary conditions (with $h=1$). As the linked page indicates, its eigenvalues are given by $$ \lambda_j = -4 \sin^2\left(\frac {\pi j}{2N} \right), \quad j = 0,1,\dots,N-1. $$ The matrix that you are considering is given by $M = A + (N-1)I$, where $I$ denotes the identity matrix. It follows that the eigenvalues of $n$ are given by $N - 1 + \lambda_j$ (with $\lambda_j$ as above), and the trace of $M^d$ is given by $$ \operatorname{Tr}(M^d) = \sum_{j=0}^{n-1} (N - 1 + \lambda_j)^d = \sum_{j=0}^{N-1} \left(N - 1 - 4\sin^2\left(\frac{\pi j}{2 N} \right)\right)^d. $$ We can rewrite this by noting that $$ \sin^2 \left(\frac{\pi j}{2N} \right) = \frac 12 \left[1 - \cos\left( \frac{\pi j}{N}\right) \right] \implies\\ N - 1 - 4\sin^2\left(\frac{\pi j}{2 N}\right) = N - 1 - \left[2 - 2\cos\left(\frac{\pi j}{N}\right)\right] = N - 3 + 2 \cos\left( \frac{\pi j}{N}\right). $$
For the case of $N = 4$, this gives us $$ \operatorname{Tr}(M^d) = \\ \left(1 + 2\cos(0)\right)^d + \left(1 + 2\cos(\pi/4)\right)^d + \left(1 + 2\cos(\pi/2)\right)^d + \left(1 + 2\cos(3 \pi /4)\right)^d = \\ 3^d + (1 + \sqrt{2})^d + 1^d + (1 - \sqrt{2})^d, $$ confirming your result. We could obtain a similar closed form for $N = 5$ by noting that $$ \cos(\pi/5) = \frac{1 + \sqrt{5}}{4}, \quad \cos(2\pi/5) = \frac{-1 + \sqrt{5}}{4}, \\ \cos(3 \pi /5) = \frac{1 - \sqrt{5}}{4},\quad \cos(4 \pi /5) = \frac{-1 - \sqrt{5}}{4}. $$
Here is an attempt to rewrite the sum. For convenience, take $p = N-3$. We have \begin{align} \sum_{j=0}^{N-1}(p + 2 \cos(\pi j/N))^d &= \sum_{j=0}^{N-1} \sum_{k=0}^d \binom dk 2^k p^{d-k}\cos^k(\pi j/N) \\ & = \sum_{k=0}^d \binom dk 2^k p^{d-k} \sum_{j=0}^{N-1}\cos^k(\pi j/N) \end{align}
I suspect that there is a relatively quick way to compute $\sum_{j=0}^{N-1}\cos^k(\pi j/N)$. In particular, write $\cos(\pi j/N) = \frac 12 (\omega^j + \omega^{-j})$, where $\omega = \exp(\pi i/N)$. From there, we have $$ \begin{align} \sum_{j=0}^{N-1}\cos^k(\pi j/N) &= \sum_{j=0}^{N-1} 2^{-k}(\omega^{j} + \omega^{-j})^k \\ & = 2^{-k} \sum_{j=0}^{N-1} \sum_{\ell = 0}^k \binom k{\ell}\omega^{(k - 2\ell)j} \\ & = 2^{-k} \sum_{\ell = 0}^k \binom{k}{\ell} \sum_{j=0}^{N-1}\omega^{(k - 2\ell)j} \\ & = 2^{-k} \sum_{\ell = 0}^k \binom{k}{\ell} \cdot \begin{cases} \frac {1 - \omega^{(k - 2\ell)N}}{1 - \omega^{(k - 2\ell)}} & 2\ell \neq k\\ N & 2\ell = k \end{cases} \\ & = 2^{-k} \sum_{\ell = 0}^k \binom{k}{\ell} \cdot \begin{cases} \frac {1 - (-1)^{(k - 2\ell)}}{1 - \omega^{(k - 2\ell)}} & 2\ell \neq k\\ N & 2\ell = k \end{cases} \end{align} $$ From there, there are two cases to consider. If $k$ is odd, then $k - 2\ell$ is odd for all $\ell$, which means that the above sum becomes \begin{align} \sum_{j=0}^{N-1}\cos^k(\pi j/N) &= 2^{-k} \sum_{\ell = 0}^k \binom{k}{\ell} \frac {1 - (-1)}{1 - \omega^{k - 2\ell}} \\ & = 2^{-k} \sum_{\ell = 0}^k \binom{k}{\ell} \frac 2{1 - \omega^{k - 2\ell}}. \end{align} Grouping the $\ell = m$ terms with the $\ell = k + 1 - m$ terms in the above some yields \begin{align} \sum_{j=0}^{N-1}\cos^k(\pi j/N) &= 2^{-k} \sum_{m = 0}^{(k-1)/2} \binom{k}{m} \left[\frac 2{1 - \omega^{k - 2m}} + \frac 2{1 - \omega^{-(k - 2m)}}\right] \\ & = 2^{-k} \sum_{m = 0}^{(k-1)/2} \binom{k}{m} 4\operatorname{Re}\left[\frac 1{1 - \omega^{k - 2m}}\right] \\ & = 2^{-k} \sum_{m = 0}^{(k-1)/2} \binom{k}{m} 4\operatorname{Re}\left[\frac {1 - \omega^{-(k - 2m)}}{|1 - \omega^{k - 2m}|^2}\right] \\ & = 2^{-k} \sum_{m = 0}^{(k-1)/2} \binom{k}{m} 4\operatorname{Re}\left[\frac {(1 - \cos((k - 2m)\pi/N)) + i\sin((k - 2m)\pi/N) }{(1 - \cos((k - 2m)\pi/N))^2 + \sin^2((k - 2m)\pi/N)}\right] \\ & = 2^{-k} \sum_{m = 0}^{(k-1)/2} \binom{k}{m} 4 \frac{1 - \cos((k - 2m)\pi/N)}{2(1 - \cos((k - 2m)\pi/N))} \\ & = 2^{-k} \sum_{m = 0}^{(k-1)/2} 2\binom{k}{m} \\ & = 2^{-k} \sum_{m = 0}^{k} \binom{k}{m} = 1 \end{align}
If $k$ is even, then the numerator $1 - (-1)^{k - 2\ell}$ will be zero for all $\ell \neq k/2$, which means that the sum becomes \begin{align} \sum_{j=0}^{N-1}\cos^k(\pi j/N) &= 2^{-k} \binom{k}{k/2} N. \end{align} That is, we have $$ \sum_{j=0}^{N-1}\cos^k(\pi j/N) = \phi(k) := \begin{cases} 2^{-k} \binom{k}{k/2} N & k \text{ is even}\\ 1 & k \text{ is odd} \end{cases} $$
Putting that together with the earlier work, we have \begin{align} \sum_{j=0}^{N-1}(p + 2 \cos(\pi j/N))^d &= \sum_{k=0}^d \binom dk 2^k p^{d-k} \sum_{j=0}^{N-1}\cos^k(\pi j/N) \\ & = \sum_{k=0}^d \binom dk 2^k p^{d-k} \phi(k) \\ & = N\left[\sum_{k \leq d, \ k \text{ even}}\binom dk \binom{k}{k/2} p^{d-k}\right] + \sum_{k \leq d, \ k \text{ odd}} \binom dk 2^k p^{d-k} \end{align} The second summation could be simplified a bit further still: \begin{align} \sum_{k \leq d, \ k \text{ odd}} \binom dk 2^k p^{d-k} &= \frac 12 \sum_{k = 0}^d \binom dk 2^k(1 - (-1)^k) p^{d-k} \\ & = \frac 12 \sum_{k = 0}^d \binom dk 2^k p^{d-k} - \frac 12 \sum_{k = 0}^d \binom dk (-2)^k p^{d-k} \\ & = \frac{(p+2)^d - (p-2)^d}{2}. \end{align}
The final formula, including corrections described below:
\begin{align} \sum_{j=0}^{N-1}(p + 2 \cos(\pi j/N))^d &= N\left[\sum_{k \leq d, \ k \text{ even}}\binom dk p^{d-k}\sum_{m = -\lfloor k/(2N)\rfloor}^{\lfloor k/(2N)\rfloor} \binom{k}{\frac k2 - Nm}\right] \\ & \qquad + \frac{(p+2)^d - (p-2)^d}{2} \end{align}
Here's a quick Python script that verifies that the (final, corrected) formula works.
The work for $\phi(k)$ for the case where $k$ is even is wrong. We instead have \begin{align} \phi(k) & = 2^{-k} \sum_{\ell = 0}^k \binom{k}{\ell} \cdot \begin{cases} \frac {1 - (-1)^{(k - 2\ell)}}{1 - \omega^{(k - 2\ell)}} & 2N \nmid (k - 2\ell)\\ N & \text{otherwise.} \end{cases} \end{align} Note that $k - 2\ell = 2Nm$ (for some integer $m$) implies that $\ell = \frac k2 - Nm$. With that, we have \begin{align} \phi(k) & = 2^{-k} \sum_{\ell = 0}^k \binom{k}{\ell} \cdot \begin{cases} \frac {1 - (-1)^{(k - 2\ell)}}{1 - \omega^{(k - 2\ell)}} & 2N \nmid (k - 2\ell)\\ N & \text{otherwise} \end{cases} \\ & = 2^{-k} N\sum_{m = -\lfloor k/(2N)\rfloor}^{\lfloor k/(2N)\rfloor} \binom{k}{\frac k2 - Nm} \end{align}