Tough Recurrence Relation

Let $b = cd$ and we will assume $b \ne 0$.

Let $f_n = f(n)$, the recurrence relation can be rewritten as

$$f_n - f_{n-2} = f_{n-1}\times \begin{cases} c, & n \text{ even }\\ -\frac{2d}{n+1},& n\text{ odd } \end{cases} \quad\iff\quad \begin{cases} f_{2m} - f_{2m-2} &= c f_{2m-1}\\ f_{2m+1} - f_{2m-1} &= -\frac{d f_{2m}}{m+1} \end{cases} $$ Let $g_m = f_{2m}$, it is clear once we figure out what $g_m$ are, we can compute $f_{2m-1}$ from $g_m$ and $g_{m-1}$.
It is clear $g_{-1} = g_0 = 1$. Furthermore, for any $m \ge 0$, we have

$$g_{m+1} - 2 g_m + g_{m-1} = (f_{2m+2}-f_{2m}) - (f_{2m} - f_{2m-2}) = c(f_{2m+1} - f_{2m-1}) = -\frac{b g_m}{m+1}$$

Let $g(z) = \sum\limits_{m=0}^\infty g_m z^m$ be the OGF. Multiply both sides of above equation by $z^{m+1}$ and start to sum from $m = 0$, we get

$$(g(z)-1) - 2zg(z) + (z + z^2g(z)) = (1-z)^2 g(z) + (z-1) = -b\int_0^z g(t)dt $$ Differentiate $z$ on both sides, we get $$\begin{align} & (1-z)^2 g'(z) - 2(1-z) g(z) + 1 = - b g(z)\\ \iff & g'(z) + \left[ -\frac{2}{1-z} + \frac{b}{(1-z)^2}\right] g = - \frac{1}{(1-z)^2}\\ \iff & \frac{d}{dz}\left[ (1-z)^2 e^{bz/(1-z)} g(z) \right] = - e^{bz/(1-z)}\\ \implies & g(z) = \frac{1}{(1-z)^2}e^{-bz/(1-z)}\left[ 1 - \int_0^z e^{bt/(1-t)} dt \right]\tag{*1} \end{align} $$ With some algebra, we can express $g(z)$ in terms of exponential integrals. $$\frac{1}{1-z} - \frac{b}{(1-z)^2} e^{-b/(1-z)}\left[{\mathrm Ei}\left(\frac{b}{1-z}\right) - {\mathrm Ei}(b)\right]$$ Unluckily, I cannot figure out how to use this to extract a nice expression of $g_n$. It turns out there is another ugly way.

Recall the definition of Laguerre polynomials $L_n(x)$ in terms of its OGF: $$\sum_{m=0}^\infty t^m L_m(x) = \frac{1}{1-t} e^{-xt/(1-t)}$$ We find

$$1 - \int_0^z e^{bt/(1-t)} dt = 1 - \int_0^1 \sum_{m=0}^\infty t^m(1-t) L_m(-b) dt = 1 - z - \sum_{m=0}^\infty \frac{L_m(-b) - L_{m-1}(b)}{m+1} z^{m+2} $$ As a result,

$$g(z) = \frac{1}{1-z}\left(\sum_{n=0}^\infty L_n(b) z^n\right)\left((1-z) - z^2 \sum_{m=0}^\infty \frac{L_m(-b) - L_{m-1}(b)}{m+1} z^m \right)$$

Compare coefficients of both sides, we find $g_n$ is a polynomial in $b$ of degree $n$.

$$f_{2n} = g_n = L_n(b) - \underbrace{\sum_{q=0}^{n-2}\frac{L_{q+1}(-b)-L_{q}(-b)}{q+2}\left(\sum_{p=0}^{n-2-q} L_p(b)\right)}_{= 0 \text{ when } n < 2}$$

Using this formula, the first few $g_{n}$ are given below $$ \begin{array}{rcl} g_0 &=& 1\\ g_1 &=& 1-b\\ 2!g_2 &=& 2-5b+b^2\\ 3!g_3 &=& 6-26b+11b^2-b^3\\ 4!g_4 &=& 24-154b+102b^2-19b^3+b^4\\ 5!g_5 &=& 120-1044b+954b^2-272b^3+29b^4-b^5 \end{array} $$ They match what one get if we solve the recurrence relations by brute force.