Incremental redundancy with ACK/NACK feedback produces a variable-length stop-feedback (VLSF) code constrained to have $m$ decoding times, with an ACK/NACK feedback to the transmitter at each decoding time. This paper focuses on the numerical evaluation of the maximal achievable rate of random VLSF codes as a function of $m$ for the binary-input additive white Gaussian noise channel, binary symmetric channel, and binary erasure channel (BEC). Leveraging Edgeworth and Petrov expansions, we develop tight approximations to the tail probability of length-$n$ cumulative information density that are accurate for any blocklength $n$. We reduce Yavas et al.'s non-asymptotic achievability bound on VLSF codes with $m$ decoding times to an integer program of minimizing the upper bound on the average blocklength subject to the average error probability, minimum gap, and integer constraints. We develop two distinct methods to solve this program. Numerical evaluations show that Polyanskiy's achievability bound for VLSF codes, which assumes $m = \infty$, can be approached with a small $m$ for all three channels. For BEC, we consider systematic transmission followed by random linear fountain coding. This allows us to obtain a new achievability bound stronger than a previous bound and new VLSF codes whose rate further outperforms Polyanskiy's bound.