A typical birth-death process is defined such as the probability of going from any state $j$ to any state $i$ is given by:
$$ p_{ij}= \begin{cases} b_i & \text {if $j = i+1$} \\ d_i & \text {if $j = i-1$} \\ 1 - b_i - d_i & \text {if $j = i$} \\ 0 & \text {otherwise}\\ \end{cases} $$
Therefore the transition matrix is:
$$ M= \left(\begin{matrix} 1-b_o & d_1 & 0 & … & 0 \\ b_0 & 1 - b_1 - d_1 & d_2 & … & 0 \\ 0 & b_1 & 1 - b_2 - d_2 & … & 0 \\ 0 & 0 & 0 & … & 1-d_m\\ \end{matrix}\right) $$
Then, one can find the stationary distribution by seeking the vector $\overrightarrow \pi$ that solves:
$$(I-M)\overrightarrow \pi = \overrightarrow 0$$
,where $\overrightarrow 0$ is a vector of zeros and $I$ is the identity matrix.
Question
If we do not respect the Markov property and every $b_i$ and $d_i$ are functions of $t$ ($b(i,t)$ and $d(i,t)$). Then, how can we calculate the stationary distribution of this pseudo-birth-death model?
In order to make sure there is a stationary distribution, I am suggesting to modify the above matrix so that when the system enters in one of the two extreme cases, it stays there forever.
$$ M= \left(\begin{matrix} 1 & d_1 & 0 & … & 0 \\ 0 & 1 - b_1 - d_1 & d_2 & … & 0 \\ 0 & b_1 & 1 - b_2 - d_2 & … & 0 \\ 0 & 0 & 0 & … & 1\\ \end{matrix}\right) $$
As an example where $b_i$ and $d_i$ are actually functions of $t$, one can consider:
$b_{it} = C_b^{\space t}$ (or $b(i,t) = C_b^{\space t}$ if you prefer) and $d_{it} = C_d^{\space t}$, where $C_b$ and $C_d$ are constant terms.
(The recurrence relation would be $b_{it} = b_{it-1} \cdot C_b$ or $b(i,t)=b(i,t-1) \cdot C_b$ if you prefer)