Updated 12 Dec 2019
PDF document and MATLAB script that demonstrates how to solve the Earth orbit, J2-perturbed form of Lambert’s problem. Lambert’s problem is concerned with the determination of an orbit that passes between two positions within a specified time-of-flight. This classic astrodynamic problem is also known as the orbital two-point boundary value problem (TPBVP) or the flyby and rendezvous problems.
David Eagle (2020). Gravity-perturbed Earth Orbit Lambert Problem - OTB (https://www.mathworks.com/matlabcentral/fileexchange/73662-gravity-perturbed-earth-orbit-lambert-problem-otb), MATLAB Central File Exchange. Retrieved .
Machine epsilon is about 2.23e-16. As you discovered Lamont, some cases may require more iterations.
I think I figured it out, the "tol" check in the continued fraction needs to be replaced by ~1e-15 and then niter needs to be allowed to hit >= 22
I think I may have found a bug in the stm2.m routine in this download (along with several other downloads that I spot checked):
[rf, vf, stm] = stm2(1, 18.96, [ 1 0 0 ]', [ -1.3673 -2.5865 -0.6116 ]')
-8.3039e+14 -9.9561e+14 -2.3542e+14
-1.6594 -1.9896 -0.47046
u0,u1,u2,u3 seem to be growing up to about 1/MAXEPS
This nearby point works fine:
[rf, vf, stm] = stm2(1, 18.96, [ 0.9995 0.00018849 0.00021848 ]', [ -1.3673 -2.5865 -0.6116 ]')
-30.729 -38.385 -9.0782
-1.6644 -1.9949 -0.47183
The elements of the bad point are:
oev = eci2orb1(1, [ 1 0 0 ]', [ -1.3673 -2.5865 -0.6116 ]')
-0.14423 7.0696 2.9094 3.6815 3.1416 5.7433