Probability of second incidence on the boundary for time t<m in Brownian motion

Question

Assume that we have a particle moving in 1D starting at $x=0$, (or call it a Martingale or Brownian Motion) where its movement is bounded at $x=-c$ and $x=c$. It moves 1 either to the left or to the right at each time slot with equal probability (no drift).

The probability that the first incidence will occur before time slot $m$ with the boundary is given as

$$ \mathrm{Pr}\{T< m\}=1-\sum_{i=0}^\infty(-1)^{i+1}\left[\Phi\left(\frac{(2i+1)c}{\sqrt{m}}\right)-\Phi\left(\frac{(2i-1)c}{\sqrt{m}}\right)\right], $$ where $\Phi()$ is the CDF of the standard normal distribution. I reached this result following a one sided boundary problem and I used reflection principle. This result is an approximation as it uses Stirling approximation to move from binomial to normal distribution, but it is as accurate as the Stirling approximation of the form

$$ n!=\sqrt{2\pi{n}}\left(\frac{n}{e}\right)^n $$

Then I started to get curious about the probability of the second incidence within time slot $m$.

One clarification, I look for hitting any boundary twice, i.e., it can be the same boundary. It is not the same as hitting opposite boundaries is within time slot $m$, which can be easily solved using $c\rightarrow 2c$. Also, I am not asking for exactly two incidences, I am asking for at least two incidences within time slot $m$.

I don't know if these will come handy, but I share the codes I used for the first incidence below.

#Performs the summation in the above equation upto ss
from scipy.stats import norm
def above_formula(c,m,ss):
    aa=c/(m)**0.5
    P1=lambda s: ((-1)**(s+1))*(norm.cdf((2*s+1)*aa)-norm.cdf((2*s-1)*aa))
    res=1
    for i in range(ss):
        res+=P1(i)+P1(-i)
    return res-P1(0)

#Simulates hitting of N particles for m time slots. 
#Returns number of particles that hasn't hit the boundary at time m
def hitting_simulation(c,m,N):
    hit=0
    for j in range(N):
        martingale=[0]
        for i in range(m):
            martingale.append(martingale[-1]+random.choice([-1,1]))
            if martingale[-1]==c or martingale[-1]==-c:
                hit+=1
                break
    return N-hit

PS: I am not the best coder out there, so please have mercy on my code.

PPS: I know that the solution is of the form $$ \sum_{i=c}^m \mathrm{Pr}\{T=i\}\mathrm{Pr^\prime}\{T<m-i\}, $$

where $\mathrm{Pr^\prime}$ is the first incidence probability of a particle starting either at $c-1$ or at $-c+1$, due to the fact that the particle has to be at either one of these locations after the first incidence. In other words, I know that we can separate the question in two parts, the second being asymmetrical first incidence problem.

Answer 1

\it's a job for markov model, assuming you choose reasonable $c\ \text{and}\ m$ you can solve it with program although with twist: idea is that we will have "game board" with size 2, 2*c-1 and particles will travel first level until they will hit a wall, then they travel to level 2 and when they hit a wall then they are added to result probability

def second_hit_prob(c, m):
    res = 0.
    space = np.full((2, 2*c-1), 0.)
    space[0, c-1]=1
    right_index = 2*c-2
    for _ in range(m):
#         print(res)
        for j in range(2):
            if j == 0: # hit left wall 1 time
                prob_wall = space[j, 0] / 2
                prob_from_left = prob_wall
                next_lvl_prob=prob_wall
                space[j, 0] = 0
            elif j == 1: # hit left wall 2 time
                prob_wall = space[j, 0] / 2
                space[j, 0] = next_lvl_prob
                res += prob_wall
                prob_from_left = prob_wall
            for i in range(1, 2*c-2):
                new_step_prob = space[j, i] / 2
                space[j, i-1] += new_step_prob
                space[j, i] = prob_from_left
                prob_from_left = new_step_prob
            if j == 0: # hit right wall 1 time
                prob_wall = space[j, right_index] / 2
                space[j, right_index-1] += prob_wall
                next_lvl_prob=prob_wall
                space[j, right_index] = prob_from_left
            elif j == 1: # hit right wall 2 time
                prob_wall = space[j, right_index] / 2
                space[j, right_index-1] += prob_wall
                space[j, right_index] = prob_from_left
                res += prob_wall
    return res
second_hit_prob(10, 300)
```