Does Kepler's Method for Determining Time of flight between two True Anomalies break down with eccentricities approaching 1?

As noted in the comments, both methods potentially break down. Kepler's method is mathematically exact, but as formulated here it becomes ill-conditioned as the eccentricity approaches 1. The mean and eccentric anomalies become undefined for a parabolic orbit because a parabola lacks a center. Numerov's/Cowell's method approximates the propagator with a Taylor series whose accuracy is not guaranteed for eccentricities all the way up to a limiting value of 1.

To separate these issues we should put Kepler's method into a form that remains well-conditioned all the way to the parabolic limit. This proper conditioning will then guarantee the accuracy of the Kepler result, against which the Numerov/Cowell result can be compared.

In this discussion the time from periapsis is used together with the anomaly conversions given below, the latter taken from Wikipedia:

$t=sqrt{dfrac{alpha^3}{mu}}M=sqrt{dfrac{r_p^3}{(1-e)^3mu}}M tag{1}label{Eq 1}$

$M=E-esin Etag{2}label{Eq 2}$

$E=2tan^{-1}left(sqrt{dfrac{1-e}{1+e}}tandfrac{theta}{2}right)tag{3}label{Eq 3}$

In Eq. 1 $r_p$ is the periapsis distance, which unlike the semi major axis remains bounded and well-defined all the way up to (and beyond) eccentricity $1$.

When $E$ is calculated from $theta$ via Eq. 3 with $e$ approaching 1, we get $E$ proportional to $(1-e)^{1/2}$; but in order to have finite and bounded times the mean anomaly $M$ must be proportional to $(1-e)^{3/2}$. Therefore Eq. 2 which connects $M$ to $E$ is ill-conditioned, for we are inputting terms with the lower power proportionality to get a difference with the higher-power proportionality.

To get $M$ in terms of quantities having the proper proportionality for a well-conditioned operation, render

$M=E-esin E = (E-sin E)+(1-e)sin E = ((sin^{-1}s)-s)+(1-e)s$

where $s=sin E$. This equation holds for $|E|lepi/2$ corresponding to $thetale 2tan^{-1}(sqrt{(1+e)/(1-e))}$. For larger true anomalies we can use Eqs. 1-3 directly because the ill-conditioning does not arise. Hereafter we focus on the case $thetale 2tan^{-1}(sqrt{(1+e)/(1-e))}$.

From Eq. 3 and the trigonometric identity $sin(2tan^{-1}u)=2u/(1+u^2)$ we obtain

$s=dfrac{2sqrt{dfrac{1-e}{1+e}}tandfrac{theta}{2}}{1+dfrac{1-e}{1+e}tan^2dfrac{theta}{2}}tag{4}label{Eq 4}$

Now we have to tackle $sin^{-1}s-s$ which is proportional to $s^3$ whereas the terms are proportional to $s$. To eliminate this ill-conditioning, convert this transcendental function to an integral involving an algebraic one:

$displaystyle{(sin^{-1}s)-s=int_0^sleft(dfrac{1}{sqrt{1-x^2}}-1right) dx = int_0^sleft(dfrac{x^2}{sqrt{1-x^2}(sqrt{1-x^2}+1)}right) dxtag{5}label{Eq 5}}$

We have invoked the difference of squares factorization to rationalize the numerator, which gets rid on the ill-conditioned subtraction. This may be integrated numerically to obtain a well-conditioned result at high eccentricity. Since the integration is independent of the orbit, values as a function of $s$ may be stored beforehand in a table for lookup.

Put this all together and we obtain the time from periapsis for a highly eccentric orbit from a modified version of Eq. 1:

$displaystyle{t=sqrt{dfrac{r_p^3}{mu}}left(dfrac{1}{(1-e)^{3/2}}int_0^sleft(dfrac{x^2}{sqrt{1-x^2}(sqrt{1-x^2}+1)}right) dx+dfrac{s}{(1-e)^{1/2}}right)}tag{6}label{Eq 6}$

with $s$ determined from Eq. 4.