[A partial answer illustrating the first method.]
State $1$ represents losing the game, so remove it from the system and condition the remaining transition probabilities on the event that the system does not transition from state $i$ to state $1$. In practical terms, this means that you delete from your transition matrix column $1$ and every row that had a $1$ in this column, and scale each remaining row $i$ by $1/(1-p_{1i})$. For your example system, this produces the reduced matrix $$P' = \begin{bmatrix}0&1&0&0&0 \\ 0&0&0&1&0 \\ q&0&0&0&p \\ 0&0&q&0&p \\ 0&0&0&0&1 \end{bmatrix}.$$ Applying standard techinques to this matrix, we can verify that the absorption probability is $1$ for any starting state and that the expected conditional absorption times are $$\begin{bmatrix}{3+q\over1-q^2} \\ {2+q+q^2\over1-q^2} \\ {1+3q\over1-q^2} \\ {1+q+2q^2\over1-q^2}\end{bmatrix}.$$
HINT: If a walk is started at 2, for it to visit 3 a second time, it must:
- travel from 2 to 3, then
- travel from 3 to 3 again.
Since expectation is linear, the time you want is the sum of the expectations of those two events, which you can find by methods like what you outlined in your question.
EDIT to clarify based on your comments: I think you're still making this much too complicated.
First, a correction to your question: you outlined a solution to find $u_2$ but I think your system of equations had some mistakes. It should have been:
\begin{align*}
u_1 &= 1 + \color{red}{\frac 1 3 u_1} + \frac 1 3 u_2 + \frac 1 3 u_4 \\
u_2 &= 1 + \frac 1 4 u_1 + \color{red}{\frac 1 4 u_2} + \frac 1 4 u_3 + \color{red}{\frac 1 4 u_4} \\
u_3 & = 0 \\
u_4 & = 0
\end{align*}
The intuition for why these relationships is valid is that from each state, you first take a single step, then weight the expected time to go from your first-step destination to 3 by the probability of each move. You can omit the $u_3$ and $u_4$ terms if you want since they're $0$, but the missing $u_1$ and $u_2$ terms will be important.
To find the expected time for the second visit somewhere, you'll need to solve some different systems of equations. If we take your example of starting the walk at 2 and computing its second visit to 3, we break that up into two occurrences; first, the walk must go from 2 to 3 a first time (taking however many steps it needs), then it must go from 3 to 3. You referenced a concern about "many cases" in your comments, but there are no cases; this is the only way to visit 3 twice.
To find the expected number of steps to go from 2 to 3, we set up a system very much like the one you outlined above:
\begin{align*}
u_1 &= 1 + \frac 1 3 u_1 + \frac 1 3 u_2 + \frac 1 3 u_4 \\
u_2 &= 1 + \frac 1 4 u_1 + \frac 1 4 u_2 + \frac 1 4 u_3 + \frac 1 4 u_4 \\
u_3 &= 0 \\
u_4 &= 1 + \frac 1 2 u_3 + \frac 1 2 u_4
\end{align*}
and note that these values of $u_1$ are different than the ones above. Solving for $u_2$ gives the first needed component.
For returning from 3 to 3, we need another system:
\begin{align*}
u'_1 &= 1 + \frac 1 3 u'_1 + \frac 1 3 u'_2 + \frac 1 3 u'_4 \\
u'_2 &= 1 + \frac 1 4 u'_1 + \frac 1 4 u'_2 + \frac 1 4 u'_3 + \frac 1 4 u_4' \\
u'_3 &= \color{blue}{\frac 1 2 \cdot 1 + \frac 1 2 u'_4} \\
u'_4 &= 1 + \frac 1 2 u'_3 + \frac 1 2 u'_4
\end{align*}
We had to adjust the $u_3$ line here so that we guarantee that the walk takes at least one step. The intuition for the blue line is that a walk started at 3 will either immediately return to 3 (in which case the excursion required just 1 move), or it will go to state 4. Note that since we only care about $u'_3$, the $u'_1$ and $u'_2$ lines are irrelevant.
Solving the first system for $u_2$ gives the expected number of moves for a walk started at 2 to visit 3. Solving the second system for $u_3'$ gives the expected number of moves for a walk started at 3 to return to 3. Adding these together gives the expected number of steps for a walk started at 2 to make its second visit to 3. Again, there are no cases to consider, because everything is baked into the existing expected values we computed.
Best Answer
Here, I am generalizing NCH's answer to this question. Consider a Markov Chain with the state space $\Omega$. I use $A$ to denote the set of absorbing states and $A^c$ to denote the set of transient states ($A\cup A^c = \Omega $). I am interested in calculating $E(V_{ij}|B_{ib})$, where the random variable $V_{ij}$ is the number of visits to State $j \in A^c$, given that the system starts from State $i \in A$, and $B_{ib}$ denotes the event for absorption at State $b \in A$ given that the system starts from State $i \in A$. We know:
$$ \Pr(V_{ij}=k|B_{ib}) = \frac{\Pr(V_{ij}=k,B_{ib}) }{\Pr(B_{ib})}. $$ The probability $\Pr(B_{ib})$ can be calculated as shown in this Wikipedia article (Subsection Absorbing Probabilities).
Let's use $T_{ij}$ to denote the event of visiting State $j$, starting from State $i$, before any absorption (not just absorption at $b$). Then $V_{ij}=k \cap B_{ib}$ includes: one time moving from $i$ to $j$, $k-1$ time moving from $j$ to $j$, and moving from $j$ to $b$ in the end without visiting $j$. That is: $$ \Pr(V_{ij}=k,B_{ib}) = \Pr(T_{ij}) \Pr(T_{jj})^{k-1} [\Pr(B_{jb})(1-\Pr(T_{jj}))] . $$ To calculate $\Pr(T_{ij})$, I will use the result in Transient Probabilities subsection of this Wikipedia article. So: $$ \begin{align} E(V_{ij}|B_{ib}) &= \sum_{k=0}^\infty k \Pr(V_{ij}=k|B_{ib}) \\ &= \sum_{k=0}^\infty k \frac{\Pr(T_{ij}) \Pr(T_{jj})^{k-1} [\Pr(B_{jb})(1-\Pr(T_{jj}))]}{\Pr(B_{ib})} \\ &= \frac{\Pr(T_{ij}) [\Pr(B_{jb})(1-\Pr(T_{jj}))]}{\Pr(B_{ib})} \sum_{k=0}^\infty k \Pr(T_{jj})^{k-1} \\ &= \frac{\Pr(T_{ij}) [\Pr(B_{jb})(1-\Pr(T_{jj}))]}{\Pr(B_{ib}) (1-\Pr(T_{jj}))^2} \\ & = \frac{\Pr(T_{ij}) \Pr(B_{jb})}{\Pr(B_{ib}) (1-\Pr(T_{jj}))}, \forall i \ne j \in A, b\in A^c. \end{align} $$
If $i = j$: $$ E(V_{ii}|B_{ib}) = \frac{1}{1-\Pr(T_{ii})}, \forall i \in A, b\in A^c. $$