This addresses the first formulation of the question, asking for the probability of absorption at $0$ or $L$ before a given time, conditionally on absorption. See below for a solution to the second formulation of the question.
One can omit the conditioning by the event that one eventually hits one of the absorbing barriers, since this event happens with full probability. Now, an often fruitful approach of this kind of problem is to look for a quasi-stationary distribution.
In words, assume that a proportion $m(n)$ of particles is at site $n$ at time zero, for every $n$ between $1$ and $L-1$, and that all the particles start to move according to the specified dynamics. Particles disappear on hitting an absorbing barrier, otherwise they continue to move from site to site. Then our goal is to have a proportion $\alpha^tm(n)$ of particles on site $n$ at time $t$, for a suitable rate $a$. If $m$ is a probability distribution, this will yield $P_m(T\ge t)=\alpha^t$ where the subscript $m$ indicates that $P(X_0=n)=m(n)$ for every $n$ between $1$ and $L-1$.
The condition of quasi-stationarity with rate $\alpha$ reads as $m(0)=m(L)=0$ and, for every site $n$ between $1$ and $L-1$,
$
\alpha m(n)=pm(n-1)+qm(n+1).
$
The general solution of this linear system of equations is $m(n)=As^n+Bt^n$ for suitable values of $A$ and $B$, where $s$ and $t$ solve the characteristic equation
$$
\alpha x=p+qx^2.
$$
The condition that $m(0)=0$ yields $B=-A$. The condition that $m(L)=0$ yields $s^L=t^L$. This can only happen if the two roots $s$ and $t$ of the characteristic equation are conjugate complex numbers whose argument $u$ is a multiple of $\pi/L$. Hence, $s=r\mathrm{e}^{\mathrm{i}u}$ and $t=r\mathrm{e}^{-\mathrm{i}u}$ where
$$
r^2=p/q,\qquad 2r\cos(u)=\alpha/q.
$$
Then $m$ is proportional to the measure $m_\alpha$ defined, for every site $n$, by
$$
m_\alpha(n)=r^n\sin(nu).
$$
The condition that $m_\alpha(n)\ge0$ for every site $n$ implies that $Lu\le\pi$. Since $Lu$ must be a multiple of $\pi$, this forces $u=\pi/L$ and
$$
\alpha=\cos(\pi/L)\sqrt{4pq}.
$$
As said at the beginning, the above proves that, for every time $t$,
$$
\sum_{n=1}^{L-1}m_\alpha(n)P_n(T\ge t)=\alpha^t|m_\alpha|,
$$
where $|m_\alpha|$ denotes the total mass of the measure $m_\alpha$. Thus, for every site $n$ and time $t$,
$$
P_n(T\ge t)\le c_n\alpha^t,\qquad c_n=|m_\alpha|/m_\alpha(n).
$$
In the other direction, one can show that, for every site $n$, there exists a positive $b_n$ such that, for every time $t$,
$$
P_n(T\ge t)\ge b_n\alpha^t.
$$
Finally, $P(T\ge t)\to0$ geometrically fast when $t\to+\infty$ with a known rate $\alpha$ of convergence. Hence,
$$
P(T\le t)=1-\Theta(\alpha^t),\qquad \alpha=\cos(\pi/L)\sqrt{4pq}.
$$
Note If $L=2$, $\alpha=0$ and one checks that $P_1(T\ge t)=0$ for every $t\ge2$. If $L=3$, $\alpha=\sqrt{pq}$ and one checks that $P_n(T\ge 2t+1)=(pq)^t$ for $n=1$ or $n=2$ and for every $t\ge0$.
This addresses the second formulation of the question, asking for the probability of absorption at $0$ before a given time, conditionally on absorption at $0$, and uses the solution above.
For every $n$, let $h(n)=P_n(T_0<T_L)$. Then $(h(n))$ is the unique solution of the linear system of equations $h(0)=1$, $h(L)=0$, and $h(n)=ph(n+1)+qh(n-1)$ for every site $n$ between $1$ and $L-1$. Let $P^0$ denote the probability $P$ conditioned by the event that the random walk hits $0$ (and as a consequence does not hit $L$). Under $P^0$, the random walk becomes the so-called Doob $h$-transform of the original random walk, in the sense that it is still a Markov chain, but with transition probabilities
$$
P^0(X_{t+1}=n+1|X_t=n)=p^0_n=1-P^0(X_{t+1}=n-1|X_t=n),
$$
where
$$
p^0_n=ph(n+1)h(n)^{-1}.
$$
As a consequence, under $P^0$, for every path $\gamma$ which starts from $\gamma(0)$, ends at $\gamma(T)$ and avoids $L$, the probability $P^0(\gamma)$ that the random walk follows $\gamma$ conditionally on $X_0=\gamma(0)$ is related to the equivalent probability $P(\gamma)$ under $P$ by the formula
$$
P^0(\gamma)=P(\gamma)h(\gamma(T))h(\gamma(0))^{-1}.
$$
Now, for every site $n$ between $1$ and $L-1$,
$$
P^0_n(T_0\ge t)=\sum_\gamma P^0(\gamma)=\sum_\gamma P(\gamma)h(\gamma(T))h(n)^{-1},
$$
where the sum enumerates all the paths $\gamma$ of length $t$ which start from $n$ and end at a site between $0$ and $L-1$ (and, as a consequence, avoid $L$). That is,
$$
P^0_n(T_0\ge t)=\sum_{k=0}^{L-1}P_n(T\ge t,X_t=k)h(k)h(n)^{-1}.
$$
It happens that the sequence $(h(k))$ is nonincreasing hence $h(L-1)\le h(k)\le1$ for every $k$ between $0$ and $L-1$ and
$$
P_n(T\ge t)h(L-1)h(n)^{-1}\le P^0_n(T_0\ge t)\le P_n(T\ge t)h(n)^{-1}.
$$
We already know that $P_n(T\ge t)\to0$ geometrically fast when $t\to+\infty$ with rate $\alpha$, hence the same asymptotics holds for $P^0_n(T_0\ge t)$.
Note finally that, for every site $n$ between $0$ and $L$,
$
h(n)=\displaystyle\frac{\left(p/q\right)^{L-n}-1}{\left(p/q\right)^L-1},
$
if $p\ne q$, and $h(n)=1-(n/L)$ otherwise, and that
$$
P_n(T\ge t)h(L-1)\le P^0_n(T_0\ge t)\le P_n(T\ge t)h(L-1)^{-1}.
$$
First let's consider a closely related enumeration problem.
Define the range of a cycle $\{x_0, x_1, x_2, \ldots, x_{k-1}, x_k=x_0\}$ to be $\max_i \left|x_{i+1}-x_i\right|$. How many cycles on $\mathbb{Z}$ are there with $k$ edges and range $\le r$? (Two cycles are equivalent if they differ by a translation or a rotation.)
If $N_k^{(r)}$ is the number of cycles with $k$ edges and range no greater than $r$, then the probability of returning safely to $0$ by making random jumps of length $\le r$ is given by
$$
P^{(r)}=\sum_{k=2}^{\infty}\frac{kN_{k}^{(r)}}{(2r)^{k}}.
$$
The factor of $k$ is needed because the starting point ($0$) can be identified with any of the nodes in a given cycle; and the factor of $(2r)^{-k}$ is the probability of taking any specific sequence of $k$ jumps. For $r=1$, there is only one legal cycle (a jump to the right and then a jump to the left), and it has two edges, so
$$
P^{(1)}=\frac{2 N_2^{(1)}}{(2r)^{2}}=\frac{2\cdot 1}{(2\cdot 1)^2}=\frac{1}{2}.
$$
For $r=2$, there are two cycles for each $k\ge 2$ (depending on whether the first rightward jump is $+1$ or $+2$), so
$$
P^{(2)}=\sum_{k=2}^{\infty}\frac{kN_k^{(2)}}{(2r)^k}=2 \sum_{k=2}^{\infty}k 4^{-k}=-2x\frac{d}{dx}\sum_{k=2}^{\infty}x^{-k}\bigg\vert_{x=4}=-2x\frac{d}{dx}\left(\frac{1}{x(x-1)}\right)\bigg\vert_{x=4}=\frac{2(2x-1)}{x(x-1)^2}\bigg\vert_{x=4}=\frac{2\cdot 7}{4\cdot 9}=\frac{7}{18}.
$$
This recovers the two known results. For $r\ge 3$, the desired enumeration can be done using a transfer matrix method. I'll make this answer "community wiki" in case someone (maybe me) has the time and energy to fill in the details for $r=3$.
The gist of the transfer matrix approach is that we can list the "states" that a particular node in the cycle can have, figure out which state-to-state transitions are allowed moving from left to right on the number line, and express the number of ways to get from the left end of a cycle to the right end (i.e., the number of different cycles) in terms of a product of many identical local matrices. In this case the "state" consists of the set of jumps to, from, and over the node, along with the range remaining to each jump. Because we're counting cycles, each node will have exactly one incoming edge and one outgoing edge. The states for $r=2$ are simple:
A blue (red) line represents a jump to the right (left), and the number next to a line indicates its remaining range. A legal transition consists of these steps:
- Decrement all numbers on the right edge by the same value (at least $1$, but no more than the smallest number on the tile).
- Select a new tile with the same number of blue and red lines on its left edge as the current tile has on its right edge.
- Place the new tile to the right of the current tile and connect lines of the same color. Make sure that all lines that pass through the new tile maintain their (decremented) values.
Legal transitions for $r=2$ are: $A$ to $(B,C,D)$ (and $A$ to $D$ can happen two ways), $B$ to $(C,D)$, and $C$ to $(B,D)$. The corresponding transfer matrix is
$$
\left(\begin{matrix}0 & 0 & 0 & 0 \\
1 & 0 & 1 & 0 \\
1 & 1 & 0 & 0 \\
2 & 1 & 1 & 0
\end{matrix}\right).
$$
Now, the number of cycles with $k$ edges will be given by
$$
N^{(2)}_k = \left(\begin{matrix}0 & 0 & 0 & 1 \end{matrix}\right)
\left(\begin{matrix}0 & 0 & 0 & 0 \\
1 & 0 & 1 & 0 \\
1 & 1 & 0 & 0 \\
2 & 1 & 1 & 0
\end{matrix}\right)^{k-1}
\left(\begin{matrix}1 \\ 0 \\ 0 \\ 0
\end{matrix}\right)
=
\left(\begin{matrix}0 & 0 & 1 \end{matrix}\right)
\left(\begin{matrix}
0 & 1 & 0 \\
1 & 0 & 0 \\
1 & 1 & 0
\end{matrix}\right)^{k-2}
\left(\begin{matrix}
1 \\ 1 \\ 2
\end{matrix}\right),
$$
where we can eliminate the first row and column of the full transfer matrix because the initial state ($A$) never recurs (i.e., the first row is all zeroes). Note that we are counting the number of ways to go from the initial state to the final state, $D$, in $k-1$ steps. In this case, the vector on the right-hand side is an eigenvector of the (reduced) transfer matrix with eigenvalue $1$, so
$$
N^{(2)}_k=\left(\begin{matrix}0 & 0 & 1 \end{matrix}\right)
\left(\begin{matrix}
1 \\ 1 \\ 2
\end{matrix}\right) = 2
$$
for any $k \ge 2$. In general, though, we will have $N_k=\sum_i c_i \lambda_i^{k-2}$, where the $\lambda_i$ are eigenvalues of the reduced transfer matrix.
The $r=2$ case can be made even simpler by noting that the set of states is symmetric under exchange of red and blue ($A$ and $D$ are themselves symmetric, and $B$ and $C$ are symmetric images of each other), so we can treat mirror pairs as single states. The transitions are then: $A$ to $D$ (two ways), $A$ to $B/C$ (two ways), $B/C$ to $B/C$, and $B/C$ to $D$. So
$$
N^{(2)}_k = \left(\begin{matrix}0 & 0 & 1 \end{matrix}\right)
\left(\begin{matrix}0 & 0 & 0 \\
2 & 1 & 0 \\
2 & 1 & 0
\end{matrix}\right)^{k-1}
\left(\begin{matrix}1 \\ 0 \\ 0
\end{matrix}\right) = \left(\begin{matrix}0 & 1 \end{matrix}\right)
\left(\begin{matrix}
1 & 0 \\
1 & 0
\end{matrix}\right)^{k-2}
\left(\begin{matrix}
2 \\ 2
\end{matrix}\right) \\ = \left(\begin{matrix}0 & 1 \end{matrix}\right)\left(\begin{matrix}
2 \\ 2
\end{matrix}\right) = 2.
$$
This additional simplification, though hardly necessary for $r=2$, will come in handy for $r=3$.
For a less trivial application of the transfer matrix approach (but still not the full $r=3$ case), I considered the asymmetric model where each jump is uniformly drawn from $\{-2,-1,1,2,3\}$. The states in this case are:
The legal transitions are: $A$ to $(B,C,D,E,H)$ (and $A$ to $H$ can happen two ways), $B$ to $(C,D,H)$ (two way from $B$ to $H$), $C$ to $(D,H)$, $D$ to $(B,H)$, $E$ to $(F,G)$, $F$ to $G$, and $G$ to $H$. So the number of cycles with $k$ edges is given by
$$
\left(\begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}\right)\left(\begin{matrix}
0 & 0 & 1 & 0 & 0 & 0 & 0\\
1 & 0 & 0 & 0 & 0 & 0 & 0\\
1 & 1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 1 & 0 & 0\\
2 & 1 & 1 & 0 & 0 & 1 & 0\\
\end{matrix}\right)^{k-2}
\left(\begin{matrix}
1 \\ 1 \\ 1 \\ 1 \\ 0 \\ 0 \\ 2
\end{matrix}\right).
$$
Using the Jordan decomposition of the reduced transfer matrix, this is
$$
\left(\begin{matrix}1 & 0 & 0 & 0 & 1 & 1 & 1 \end{matrix}\right)\left(\begin{matrix}
0 & 1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 1 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & \nu_1 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & \nu_2 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & \nu_2^{*}
\end{matrix}\right)^{k-2}\left(\begin{matrix}
-1 \\ 0 \\ 1 \\ 1 \\ c_1 \\ c_2 \\ c_2^{*}
\end{matrix}\right) \\ =c_1 \nu_1^{k-2} + c_2 \nu_2^{k-2} + c_2^{*} (\nu_2^{*})^{k-2} - \delta_{k,2} + \delta_{k,4} + \delta_{k,5},
$$
where $\nu_1=1.32472$, $\nu_2=-0.662359 - 0.56228 i$, $c_1=2.945975$, and $c_2=0.0270124 + 0.118443 i$. Then the probability of safe return to $0$ is given by
$$
P^{(3,2)}=\sum_{k=2}^{\infty}\frac{kN_k^{(3,2)}}{5^k}=-\frac{2}{5^2}+\frac{4}{5^4}+\frac{5}{5^5}+\sum_{i} c_i \sum_{k=2}^{\infty}\frac{k \nu_i^{k-2}}{5^k} \\ = -\frac{9}{125} + \sum_i \frac{c_i(10-\nu_i)}{5(5-\nu_i)^2} = -0.072 + 0.378409 + 2\cdot 0.00289355 \\ = 0.3121961.
$$
For comparison, in a simple simulation of $10^8$ trials (running each trial for at most $1000$ jumps), $31218647$ asymmetric jumpers returned safely to the origin. This gives an estimate of $0.31219 \pm 0.00006$ for the probability; i.e., the simulation agrees perfectly with the matrix calculation above.
Finally, for $r=3$ we have the following set of states:
Several pairs of states are related by red-blue exchange symmetry (e.g., $B$ and $B'$). Using this symmetry, we can write the transfer matrix in terms of symmetric states only (e.g., $B/B'$), as described earlier for the $r=2$ case. The transitions out of unprimed states are:
- $A$ to $(B,B',C,C',D,H(\times 3))$
- $B$ to $(B',C,C',E,H(\times 2))$
- $C$ to $(B',H)$
- $D$ to $(F,F',G,G')$
- $E$ to $(F',G')$
- $F$ to $G$
- $G$ to $(C,H)$.
So
$$
N_k^{(3)} = \left(\begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}\right)\left(\begin{matrix}
1 & 1 & 0 & 0 & 0 & 0 & 0\\
2 & 0 & 0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0\\
1 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 2 & 1 & 0 & 0 & 0\\
0 & 0 & 2 & 1 & 1 & 0 & 0\\
2 & 1 & 0 & 0 & 0 & 1 & 0\\
\end{matrix}\right)^{k-2}
\left(\begin{matrix}
2 \\ 2 \\ 1 \\ 0 \\ 0 \\ 0 \\ 3
\end{matrix}\right).
$$
This sequence begins with $3,6,14,30,62,130,\ldots$. After diagonalizing the matrix, we arrive at
$$
N_k^{(3)}=\delta_{k,2}+\left(\begin{matrix}1 & 1 & 1 & 1 \end{matrix}\right)\left(\begin{matrix}
\nu_1 & 0 & 0 & 0 \\
0 & \nu_2 & 0 & 0 \\
0 & 0 & \nu_3 & 0\\
0 & 0 & 0 & \nu_3^{*}
\end{matrix}\right)^{k-2}\left(\begin{matrix}
c_1 \\ c_2 \\ c_3 \\ c_3^{*}
\end{matrix}\right)=\delta_{k,2} + \sum_i c_i \nu_i^{k-2},
$$
where
$$
\nu_1 = -0.716673 \\
\nu_2 = 2.10692 \\
\nu_3 = 0.304877 - 0.754529 i \\
c_1 = -0.188335 \\
c_2 = 3.136179 \\
c_3 = -0.4739206 - 0.3006352 i.
$$
Finally, the desired probability is
$$
P^{(3)} = \sum_{k=2}^{\infty}\frac{kN_k^{(3)}}{6^k}=\frac{1}{18}+\sum_i \frac{c_i(12-\nu_i)}{6(6-\nu_i)^2}=0.325873.
$$
Checking by simulation for $10^8$ trials yields an estimate of $0.32587 \pm 0.00006$, wholly consistent with this result.
Automating the process
One can do all of the following by making a number of formalizations.
It is advantageous to mechanize the process before proceeding to larger cases; in particular, we will represent a state as a triple $(X,Y,P)$ where $X$ is the set of remaining ranges of rightward (red) edges, $Y$ is the set of leftward (blue) edges, and $P$ is a partition of $X\cup Y$ representing which edges have been connected (noting that if $x=y$ then the edges represented by $x\in X$ and $y\in Y$ originated together and are clearly connected). We will say $(X',Y',P')$ is an $n$-translate of $(X,Y,P)$ if we might see the former after moving our "window" $n$ steps to the right. For convenience, define
$$m=\min(X\cup Y)$$
$$S-n=\{s-n:s\in S\}$$
and, because the relations are not conducive to standard notation, I will leave to the reader the problem of relating $P'$ and $P$. Then, we can have the following situations.
Birth: In the state where the node in question has two rightward edges passing from it. We require:
$$X'=X-n\cup \{r\}$$
$$Y'=Y-n\cup \{r\}$$
We will consider $r$ to be its own connectedness class in $P'$, and otherwise just shift everything.
Right-Continuation: This is the state where the node receives a rightward edges and also outputs one, we need that for some $c\in X-n$ we have:
$$X'=X-n\cup\{r\}\setminus\{c\}$$
$$Y'=Y-n$$
The partition will essentially shift and add $r$ to the equivalence class of $c$.
Left-Continuation: This is where a node receives and produces a leftward edge. Flop $X$ and $Y$ in the above definition.
Death: This is where a node is a local maximum, receiving two edges. We need that
$$X'=X-n\setminus\{c_x\}$$
$$Y'=Y-n\setminus\{c_y\}$$
where $c_x$ and $c_y$ are appropriate members of $X-n$ and $Y-n$ respectively. Neither $X'$ nor $Y'$ may contain $0$. The equivalence classes representing $c_x$ and $c_y$ will be merged. If $c_x$ and $c_y$ were the last remaining representatives of their equivalence class, then the remaining edges were not connected to them; this transition is invalid unless $X'=Y'=\{\}$.
Then, in the matrix, the number of ways to transition from one state to another is the number of $n$'s such that the latter is an $n$-translate of the former. My code works by generating every conceivable state - i.e. every pair $(X,Y,P$) where $|X|=|Y|$ and $P$ is a partition of $X\cup Y$. Some of these are inaccessible from the starting state or cannot access the ending state, but I did not write code to prune such states away. I also did not take advantage of the symmetry $(X,Y,P)\rightarrow(Y,X,P)$. Either would likely lead to big increases in efficiency, allowing the computation of more values $P^{(r)}$.
It is worthy of note that
$$\sum_{k=2}^{\infty}a\cdot \frac{kM^{k-1}}{2k}\cdot b=a\cdot\left(\sum_{k=2}^{\infty}\frac{kM^{k-1}}{2r}\right)\cdot b$$
by linearity, meaning we can first use the identity that
$$\sum_{k=2}^{\infty}\frac{kM^{k-1}}{(2r)^k}=\frac{1}{4r^2}(2I-\frac{1}{2r}M)\cdot M\cdot(\frac{1}{2r}M-I)^{-2}$$
which is easily derived by noting that $M$ acts, for these purposes, just like a real number with $|M|<2r$ as $M$ must have no eigenvalues of absolute value $2r$ or more (given that there are only $(2r)^k$ paths of length $k$, whereas there are about $k\lambda^k$ cycles of length $k$ with identified starting point where $\lambda$ is the largest eigenvalue of $M$). Then, we can take the transfer matrix, directly compute the infinite sum as above, and extract the desired coefficient. It should be noted that this will always yield a rational answer. In particular, we receive:
$$P^{(1)}=\frac{1}2$$
$$P^{(2)}=\frac{7}{18}$$
$$P^{(3)}=\frac{4368595}{13405842}\approx 0.325872$$
$$P^{(4)}=\frac{33373525946827013707062414432976}{117177232862732965149684560993569}\approx 0.284812$$
Mathematica code may be found in this revision. It may be noted that the methods may be modified easily to handle asymmetric cases or cases with non-uniform probability - essentially, instead of using our transfer matrix as a count, we can use it to measure the probability of a cycle occurring given that we start on a given point of it.
Best Answer
A standard approach for a probabilist here is to rely on generating functions. That is, fix $|s|<1$ and consider $u_x=\mathrm E_x(s^T)$, where the subscript $x$ means that the random walk starts from $x$ and where $T$ denotes the first return to $0$ of the random walk, that is, $T=\inf\{n\ge1\mid X_n=0\}$ and $(X_n)_n$ denotes the random walk starting from $X_0=x$ and making steps $+1$ with probability $p$ and $-1$ with probability $q=1-p$.
What you are asking amounts to computing $u_0$ and $u_1$, we shall explain how to compute $u_x$ for every integer $x$.
Taking into account the first step of the walk, one sees that $u_0=s(pu_1+qu_{-1})$, $u_1=s(pu_2+q)$ and $u_{-1}=s(p+qu_{-2})$. Furthermore, starting from $x=2$, to reach $0$ one must first hit $1$ and then starting from $1$ one must hit $0$. The times to hit $1$ starting from $2$ and to hit $0$ starting from $1$ to are i.i.d. hence $u_2=u_1^2$. Likewise, $u_{-2}=u_{-1}^2$ (and in fact $u_x=u_1^x$ and $u_{-x}=u_{-1}^x$ for every positive $x$).
This shows that $u_0=s(pu_1+qu_{-1})$ where $u_1=s(pu_1^2+q)$ and $u_{-1}=s(p+qu_{-1}^2)$. Solving this for $u_1$ and $u_{-1}$ yields $u_1=\dfrac{1\pm\sqrt{1-4pqs^2}}{2sp}$ and a similar formula for $u_{-1}$ but since one wants that $u_1$ and $u_{-1}$ stay bounded when $s\to0$, the $\pm$ signs are in fact $-$ signs and $$ u_0=1-\sqrt{1-4pqs^2},\qquad u_1=\frac{u_0}{2sp},\qquad u_{-1}=\frac{u_0}{2sq}. $$ The PDF of $T$ starting from $0$, $1$ and $-1$ is fully encoded in $u_0$, $u_1$ and $u_{-1}$ respectively. For example, the probability to come back at $0$ starting from $0$ and to hit $0$ starting from $1$ and $-1$ respectively, are the limits of $u_0$, $u_1$ and $u_{-1}$ when $s\to1$, that is, $$ \mathrm P_0(T<\infty)=1-\sqrt{1-4pq}=1-|1-2p|, $$ and $$ \mathrm P_1(T<\infty)=\frac1{2p}\mathrm P_0(T<\infty),\quad \mathrm P_{-1}(T<\infty)=\frac1{2q}\mathrm P_0(T<\infty). $$ Let us assume from now on that $p\ge\frac12$ (hence $p\ge q$). Then, $$ \mathrm P_0(T<\infty)=2q,\quad \mathrm P_1(T<\infty)=q/p,\quad\mathrm P_{-1}(T<\infty)=1. $$ Likewise, to get the full PDF only requires to know the expansion along powers of $t$ of the square root involved, namely, $$ \sqrt{1-4t}=1-\sum\limits_{n=1}^{+\infty}c_nt^n,\quad c_n=\frac{(2n)!}{(2n-1)(n!)^2}. $$ One gets for example $$ u_1=\frac1{2ps}\sum\limits_{n=1}^{+\infty}c_n(pq)^ns^{2n}, $$ hence, for every $n\ge1$, $$ \mathrm P_1(T=2n-1)=\frac1{2p}c_n(pq)^n. $$ Likewise, for every $n\ge1$, $$ \mathrm P_{-1}(T=2n-1)=\frac1{2q}c_n(pq)^n,$$ and finally, $$ \mathrm P_0(T=2n)=c_n(pq)^n. $$