How to compute the variance of number of coin flips to see HTH sequence using linear system of equations. Assuming fair coin flips, I know how to compute the expected number of coin flips to see HTH sequence by writing out the linear system of equations from the state transition diagram below. 

Define $E(X)$ as the expected number of steps from state X to state HTH. E(0) means expected number of steps to get from state 0 to state HTH and E(H) means expected number of steps to get from state H to state HTH, etc. I'd like out the equations with the transition probability and expected value of the neighbor state plus one for the immediate step. The full linear system of equations is as follows.
\begin{align}
E(0) &= \frac{1}{2} (E(H) + 1) + \frac{1}{2} (E(0) + 1)\\
E(H) &= \frac{1}{2} (E(H) + 1) + \frac{1}{2} (E(HT) + 1)\\
E(HT) &= \frac{1}{2} (E(HTH) + 1) + \frac{1}{2} (E(0) + 1)\\
E(HTH) &= 0
\end{align}
We can solve for $E(0)$ since we have 4 equations and 4 unknowns. It turns out $E(0)=10$
The question is how do we compute the variance $V(0)$. I understand that we can build the transition matrix and use the one line matrix formula from Wikipedia. However, I'm interested in intuitive system of equations approach (without any of Markov chain jargons like fundamental matrix, ergodic, transient, etc) just like the one I wrote above to compute $E(0)$. My hunch is we need systems of equations defined with $E(X^2)$ variables but it's not clear to me how to approach that.
 A: Define $F(X)$ as the expected (#steps$^2$) from state $X$ to state $HTH$. Then adding $1$ step in a recursive definition of $F$ can be done as follows:
\begin{eqnarray*}
F(0) &=& \dfrac{1}{2}\sum_k{(k+1)^2P(\text{#Steps from $H$ to $HTH$}=k)}
 + \\
&&\qquad\dfrac{1}{2}\sum_k{(k+1)^2P(\text{#Steps from $0$ to $HTH$}=k)} \\
&& \\
&=& \dfrac{1}{2}\sum_k{(k^2+2k+1)P(\text{#Steps from $H$ to $HTH$}=k)}
 + \\
&&\qquad\dfrac{1}{2}\sum_k{(k^2+2k+1)P(\text{#Steps from $0$ to $HTH$}=k)} \\
&& \\
&=& \dfrac{1}{2}\left(F(H) + 2E(H) + 1\right) + \dfrac{1}{2}\left(F(0) + 2E(0) + 1\right).
\end{eqnarray*}
The other equations are derived similarly, so we have the system of equations:
\begin{eqnarray*}
F(0) &=& \dfrac{1}{2}\left(F(H) + 2E(H) + 1\right) + \dfrac{1}{2}\left(F(0) + 2E(0) + 1\right) \\
F(H) &=& \dfrac{1}{2}\left(F(H) + 2E(H) + 1\right) + \dfrac{1}{2}\left(F(HT) + 2E(HT) + 1\right) \\
F(HT) &=& \dfrac{1}{2}\left(F(HTH) + 2E(HTH) + 1\right) + \dfrac{1}{2}\left(F(0) + 2E(0) + 1\right) \\
F(HTH) &=& 0.
\end{eqnarray*}
The system of equations in the question provides:
$$E(0)=10,\quad E(H)=8,\quad E(HT)=6.$$
So solving the system of equations gives $F(0)=114$. So the variance = $114-10^2 = 14$.
A: One approach might be to find the probability generating function, which turns out (based on  OEIS A005314) to be $$G(t)=\frac{{t}^{3}}{8-8t+2t^2-t^3}.$$
You can then calculate the mean as $G'(1^-)$, which is indeed $10$, and the variance as $G''(1^-)+G'(1^-) - (G'(1^-))^2$, which is $58$.
