Is there something like a normal distribution model for discrete probability?
Your ultimate goal is not clear. Perhaps I can flounder around and make some useful comments.
For appropriate choices of $n$ and $\theta,$ the distribution $Binom(n, \theta)$ is approximately normal, especially if $n$ is large and $\theta$ is not too far from 1/2. The mean is $\mu = n\theta$ and the variance is $\sigma^2 = n\theta(1-\theta).$
Also, for large enough $\lambda,$ the distribution $Pois(\lambda)$ is nearly normal. The mean and variance are $\mu = \lambda$ and $\sigma^2 = \lambda.$ However, the Poisson model may have less flexibility in matching what you want.
Of course, to find the probability that a random variable taking integer values lies in an interval $(a, b]$ you will add probabilities for integer values in that interval, rather than evaluating an integral.
For example, if $X \sim Binom(n = 100, \theta = 1/2),$ you have $\mu = 50$ and $\sigma = 5.$ Perhaps you want
$$P(48 < X \le 52) = P(X = 49) + P(X = 50) + P(X = 51) + P(X = 52)\\ = P(X \le 52) -P(X \le 48) = F_X(52) - F_X(48) = 0.3091736,$$ where $F_X(\cdot)$ is the CDF of $X.$
If there are many integers in the desired interval, computation by hand
can be tedious. In R statistical software dbinom
denotes a binomial PDF
and pbinom
a binomial CDF.
The probability above could be evaluated in R as shown below. [The last value is a normal approximation (with continuity correction), which is often accurate to a couple of decimal places.]
sum(dbinom(49:52, 100, .5)) # adding terms of the PDF
## 0.3091736
diff(pbinom(c(48,52), 100, .5)) # subtracting two CDF values
## 0.3091736
diff(pnorm(c(48.5,52.5), 50, 5)) # normal approximation
## 0.3093739
The figure below shows several values of the PDF of $Binom(100, .5),$ emphasizes the four probabilities required (heights of thick blue bars), and shows the approximating normal density curve. The normal approximation is the area beneath the curve between the vertical green lines.
You can construct a discrete distribution called the Skellam distribution (https://en.wikipedia.org/wiki/Skellam_distribution) from the difference of two Poisson distributions.
If $X \sim Pois(\lambda_1)$ and $Y \sim Pois(\lambda_2)$, then the distribution for $Z = X-Y$ fits the Skellam distribution with mean $\lambda_1-\lambda_2$ and variance $\lambda_1+\lambda_2$. The distribution is approximately normal with the PDF given by $Pr\{Z=n\} = e^{-(\lambda_1+\lambda_2)}\left(\sqrt{\frac{\lambda_1}{\lambda_2}}\right)^n I_n\left(2\sqrt{\lambda_1\lambda_2}\right)$.