# Generalization of Tauchen 1986 approach to a case of time-varying volatility

My question is about generalization of Tauchen'86 approach to a case of time-varying volatility.

Say, I have a process $$z_{t+1}=\rho z_t+\sigma_t \varepsilon_{t+1}$$ where $\varepsilon\sim \mathcal{N}(0,1)$ and $\sigma_t\in\{\sigma_L,\sigma_H\}$ with $\sigma_L<\sigma_H$.

I want to get the nodes and a corresponding transition matrix. Does anyone know how to implement that in Matlab?