Distribution in Heston
$$ dV_t=-k(V_t-1)dt+ \epsilon\sqrt{V_t}dW_t $$ $ W_t $ is wiener process and the rest is just some parameters.
For $ T_{i+1}>T_{i} $ how do I find the expectation and variance of $ V_{T_{i+1}} $ conitional to $ V_{T_i} $ ?
As per @Quantuple, we may consider a more general Cox-Ingersoll-Ross model
$$ {\rm d}r_t=a\left(b-r_t\right){\rm d}t+\sigma\sqrt{r_t}{\rm d}w_t. $$ Integrating this equation on $ t\in\left(u,v\right)\subseteq\mathbb{R}^+ $ yields $$ \begin{align} r_v-r_u=\int_u^v{\rm d}r_t&=\int_u^v\left[a\left(b-r_t\right){\rm d}t+\sigma\sqrt{r_t}{\rm d}w_t\right]\ &=ab\left(v-u\right)-a\int_u^vr_t{\rm d}t+\sigma\int_u^v\sqrt{r_t}{\rm d}w_t. \end{align} $$ Acting $ \mathbb{E}\left(\cdot|r_u\right) $ on both sides of this equality gives $$ \begin{align} \mathbb{E}\left(r_v|r_u\right)-r_u&=\mathbb{E}\left(r_v-r_u|r_u\right)=\mathbb{E}\left[ab\left(v-u\right)-a\int_u^vr_t{\rm d}t+\sigma\int_u^v\sqrt{r_t}{\rm d}w_t\Bigg|r_u\right]\ &=ab\left(v-u\right)-a\mathbb{E}\left(\int_u^vr_t{\rm d}t\Bigg|r_u\right)+\sigma\mathbb{E}\left(\int_u^v\sqrt{r_t}{\rm d}w_t\Bigg|r_u\right)\ &=ab\left(v-u\right)-a\mathbb{E}\left(\int_u^vr_t{\rm d}t\Bigg|r_u\right)\ &=ab\left(v-u\right)-a\int_u^v\mathbb{E}\left(r_t|r_u\right){\rm d}t. \end{align} $$ Denote $ f(t)=\mathbb{E}\left(r_t|r_u\right) $ for $ t\ge u $ , and the last equation is equivalent to, for all $ v\ge u $ , $$ f(v)=r_u+ab\left(v-u\right)-a\int_u^vf(t){\rm d}t. $$ Regard this as an integral equation with respect to $ v $ , and its respective differential equation with initial condition reads $$ \begin{align} f’(v)&=ab-af(v),\ f(u)&=r_u. \end{align} $$ This system immediately leads to
$$ \mathbb{E}\left(r_v|r_u\right)=f(v)=r_ue^{-a\left(v-u\right)}+b\left(1-e^{-a\left(v-u\right)}\right). $$
Based on the result above, we may figure out $ \mathbb{E}\left(r_v^2|r_u\right) $ as follows. Thanks to Ito’s formula, the Cox-Ingersoll-Ross model gives
$$ \begin{align} {\rm d}\left(r_t^2\right)=2r_t{\rm d}r_t+{\rm d}\left<r\right>_t&=2ar_t\left(b-r_t\right){\rm d}t+\sigma r_t\sqrt{r_t}{\rm d}w_t+\sigma^2r_t{\rm d}t\ &=\left[\left(2ab+\sigma^2\right)r_t-2ar_t^2\right]{\rm d}t+\sigma r_t\sqrt{r_t}{\rm d}w_t. \end{align} $$ Again, Integrating this equation on $ t\in\left(u,v\right)\subseteq\mathbb{R}^+ $ yields $$ r_v^2-r_u^2=\left(2ab+\sigma^2\right)\int_u^vr_t{\rm d}t-2a\int_u^vr_t^2{\rm d}t+\sigma\int_u^vr_t\sqrt{r_t}{\rm d}w_t. $$ Acting $ \mathbb{E}\left(\cdot|r_u\right) $ on this equality yields $$ g(v)-r_u^2=\left(2ab+\sigma^2\right)\int_u^vf(t){\rm d}t-2a\int_u^vg(t){\rm d}t, $$ where $ f $ agrees with the notation from above, while $ g(v)=\mathbb{E}\left(r_t^2|r_u\right) $ . Note that its respective differential equation reads $$ \begin{align} g’(v)&=\left(2ab+\sigma^2\right)f(v)-2ag(v),\ g(u)&=r_u^2. \end{align} $$ Thanks to this system, we could eventually figure out
$$ \begin{align} \text{Var}\left(r_v^2|r_u\right)&=\mathbb{E}\left(r_v^2|r_u\right)-\mathbb{E}^2\left(r_v|r_u\right)=g(v)-f^2(v)\ &=\frac{\sigma^2}{a}\left[r_ue^{-a\left(v-u\right)}\left(1-e^{-a\left(v-u\right)}\right)+\frac{b}{2}\left(1-e^{-a\left(v-u\right)}\right)^2\right]. \end{align} $$
The answer to your original question could be obtained by taking $ a=k $ , $ b=1 $ , $ \sigma=\epsilon $ , $ u=T_i $ and $ v=T_{I+1} $ .