I am working with a scalar stochastic process $X_t$ which satisfies some stochastic differential equation of the form $$ dX_t = \mu(X_t,t)dt+\sigma(X_t,t)dB_t $$ and for which I managed to compute a probability distribution $p(x,t)$ satisfying the corresponding Kolgomorov's backward equation: $$ -\partial_t p(x,t) = \mu(x,t)\partial_x p(x,t)+\frac{1}{2}\sigma(x,t)^2\partial_x^2p(x,t) $$ From there, given $t=t_0$ I know that I can compute the probability of $X_{t_0}$ to lie in a certain region using $p(x,t_0)$ as a probability distribution density. However, I am confused about what happens if I want to consider other options for $t$ other than a single value.
For example, I wonder if I can start answering questions such as computing the probability of $X_t$ being identically 0 for all $t\geq 0$, or (which I am more interested in) the probability of $X_t\in[x_0,x_1],\forall t\in[t_0,t_1]$ for given $x_0,x_1,t_0,t_1\in\mathbb{R}$. I cannot finish to understand how the information of $p(x,t)$ or the stochastic differential equation actually captures those probabilities.
EDIT: In case it helps, consider constant coefficients in the SDE and initial condition $X_0$ to be gaussian, such that $X_t$ is a guassian process. Even in this case, I am not sure how to approach the aforementioned questions.
Knowing the probability distribution density is not enough to answer the questions you are interested in. Instead, you may write an equation for some quantity of interest.
For example, let us consider the probability that $X_t \in [a,b]$ for all $t\in [0,T]$. To write an equation, denote $\tau_{[a,b]} = \inf\{s: X_s\notin [a,b]\}\wedge T$ and define $$ q(x,t) = \mathbb{P}(X_s\in [a,b]\ \forall s\in[t,T]\mid X_t = x) = \mathbb{P}\big(\tau_{[a,b]}\ge T\mid X_t = x\big)\\ = \mathbb{E}[ \mathbb{1}_{\tau_{[a,b]}\ge T}\mid X_t = x] = \mathbb{E}\big[ g({X}_{\tau_{[a,b]}})\mid X_t = x\big], $$ where $g(x) = \mathbb{1}_{(a,b)}(x)$. Despite $g$ is discontinuous, we can still write the boundary value problem for $q$: $$ \begin{cases} \frac{\partial}{\partial t} q(x,t) + \mathcal L_{t,x} q(x,t) = 0,\quad & x\in (a,b), t<T;\\ q(x,T) = 1,& x \in (a,b);\\ q(a,t) = q(b,t) = 0, & t<T, \end{cases} $$ where $\mathcal L_{t,x} = \mu(x,t)\partial_x +\frac{1}{2}\sigma(x,t)^2\partial_x^2$ is the generator of $X$.
This boundary value problem is quite hard to solve in general, especially for time-dependent $\mu$ and $\sigma$. In homogeneous case, more or less explicit solutions for some particular cases (Brownian motion with or without drift, Ornstein–Uhlenbeck process) can be found in Borodin and Salminen's "Handbook".