Errores de trayectoria del solucionador Lambert de tiempo de vuelo bajo

Tengo un solucionador de Lambert que se usa para calcular el delta-V requerido para las trayectorias de transferencia interplanetaria, pero recientemente me encontré con un problema. Funciona bien para trayectorias de tiempo de vuelo "largo" (por ejemplo, una trayectoria de 250 días de la Tierra a Marte), pero si pruebo una trayectoria de tiempo de vuelo "corto" (por ejemplo, una trayectoria de 45 días desde Earth to Mars), el solucionador parece romperse, al no converger sus cálculos. Me doy cuenta de que un TOF de 45 días es bastante poco realista, pero generalmente realizo una corrección a mitad de camino a la mitad de una trayectoria, por lo que si tengo un TOF rápido de 90 días, realizaré un MCC 45 días en el vuelo (que es cuando el solucionador se rompe).

El siguiente es el solucionador de Lambert que utiliza una formulación de variables universales (codificada en el lenguaje de Mathematica) y un solucionador de bisección iterativo simple para calcular el delta-V requerido dado un tiempo de vuelo requerido:

EDITAR: desde entonces arreglé el código y el solucionador de Lambert a continuación parece funcionar bastante bien:

G = 6.672*10^-11; (*Gravitational Constant*)
m[0] = 1.988544*10^30; (*Mass of Sun*)
TOF = (45) (86400); (*Time-of-flight in seconds*)
R[1] = {-1.1751563715176448`*^11, 8.982523733108002`*^10} (*Heliocentric position of Earth*)
R[2] = {-5.256112524631399`*^10, -2.1604439066188406`*^11} (*Heliocentric position of Mars*)
R1 = Sqrt[R[1].R[1]]
R2 = Sqrt[R[2].R[2]]
\[CapitalDelta]\[Nu] = ArcCos[R[1].R[2]/(R1 R2)] (*Change in true anomaly, in radians*)
A = Sqrt[R1 R2 (1 + Cos[\[CapitalDelta]\[Nu] ])];

iterationCount = 0;
z = 0;
zhi = 4 \[Pi]^2;
zlow = -4 \[Pi];
c[z_] := If[z > 0, (1 - Cos[Sqrt[z]])/z, 
  If[z < 0, (1 - Cosh[Sqrt[-z]])/z, 1/2]]
S[z_] := If[z > 0, (Sqrt[z] - Sin[Sqrt[z]])/Sqrt[z^3], 
  If[z < 0, (Sinh[Sqrt[-z]] - Sqrt[-z])/Sqrt[(-z)^3], 1/6]]
Y[z_] := R1 + R2 - (A (1 - S[z] z))/Sqrt[c[z]];
X[z_] := Sqrt[Y[z]/c[z]];
t[z_] := (X[z]^3 S[z] + A Sqrt[Y[z]])/Sqrt[G m[0]];

t[z] = t[z]; (*Initial value for t[z]*)

(*Iterative Bisection Solver*)
While[Norm[t[z] - \[CapitalDelta]t] > 1*10^-6 && iterationCount < 100,
 c[z_] := 
  If[z > 0, (1 - Cos[Sqrt[z]])/z, 
   If[z < 0, (1 - Cosh[Sqrt[-z]])/z, 1/2]];
 S[z_] := 
  If[z > 0, (Sqrt[z] - Sin[Sqrt[z]])/Sqrt[z^3], 
   If[z < 0, (Sinh[Sqrt[-z]] - Sqrt[-z])/Sqrt[(-z)^3], 1/6]];
 Y[z_] := R1 + R2 + (A (S[z] z - 1))/Sqrt[c[z]];
 (*Making sure Y>0*)
 While[A > 0 && Y[z] < 0,
  zlow = zlow + 0.01;
  z = (zhi + zlow)/2;
  c[z_] := 
   If[z > 0, (1 - Cos[Sqrt[z]])/z, 
    If[z < 0, (1 - Cosh[Sqrt[-z]])/z, 1/2]];
  S[z_] := 
   If[z > 0, (Sqrt[z] - Sin[Sqrt[z]])/Sqrt[z^3], 
    If[z < 0, (Sinh[Sqrt[-z]] - Sqrt[-z])/Sqrt[(-z)^3], 1/6]];
  Y[z_] := R1 + R2 + (A (S[z] z - 1))/Sqrt[c[z]];
  ];
 X[z_] := Sqrt[Y[z]/c[z]];
 t[z_] := (X[z]^3 S[z] + A Sqrt[Y[z]])/Sqrt[G m[0]];
 If[t[z] <= \[CapitalDelta]t, zlow = z, zhi = z];
 z = (zhi + zlow)/
  2; (*Re-calculating z using bisection root finding method*)
 Print[Norm[t[z] - \[CapitalDelta]t]];
 iterationCount++;];

    f = 1 - Y[z]/R1;
    g = A Sqrt[Y[z]/(G m[0])];
    gdot = 1 - Y[z]/R2;
    vLambert[1] = (R[2] - f R[1])/g; (*Required velocity at start of trajectory*)
    vLambert[2] = (gdot R[2] - R[1])/g; (*Velocity at target arrival*)

Como se puede ver, existen controles para asegurarse de que Y siempre se mantenga positivo, por lo que no estoy muy seguro de cuál podría ser el problema. Desafortunadamente, buscar en Internet artículos académicos o capítulos de libros tampoco me ha dado una solución a este problema. Cualquier ayuda sería muy apreciada, muchas gracias.