For further understanding of neural coding, stochastic variability of interspike intervals has been investigated by both experimental and theoretical neuroscientists. In stochastic neuronal models, the interspike interval corresponds to the time period during which the process imitating the membrane potential reaches a threshold for the first time from a reset depolarization. For neurons belonging to complex networks in the brain, stochastic diffusion processes are often used to approximate the time course of the membrane potential. The interspike interval is then viewed as the first passage time for the employed diffusion process. Due to a lack of analytical solution for the related first passage time problem for most diffusion neuronal models, a numerical integration method, which serves to compute first passage time moments on the basis of the Siegert recursive formula, is presented in this paper. For their neurobiological plausibility, the method here is associated with diffusion processes whose state spaces are restricted to finite intervals, but it can also be applied to other diffusion processes and in other (non-neuronal) contexts. The capability of the method is demonstrated in numerical examples and the relation between the integration step, accuracy of calculation and amount of computing time required is discussed.