diff --git a/docs/source/examples/revisiting-lamberts-problem-in-python.myst.md b/docs/source/examples/revisiting-lamberts-problem-in-python.myst.md
index d931ea226..b96a1cf6f 100644
--- a/docs/source/examples/revisiting-lamberts-problem-in-python.myst.md
+++ b/docs/source/examples/revisiting-lamberts-problem-in-python.myst.md
@@ -4,9 +4,9 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.0
+    jupytext_version: 1.16.0
 kernelspec:
-  display_name: Python 3
+  display_name: Python 3 (ipykernel)
   language: python
   name: python3
 ---
@@ -15,25 +15,25 @@ kernelspec:
 
 The Izzo algorithm to solve the Lambert problem is available in hapsira and was implemented from [this paper](https://arxiv.org/abs/1403.2705).
 
-```{code-cell}
+```{code-cell} ipython3
 from cycler import cycler
 
 from matplotlib import pyplot as plt
 import numpy as np
 
-from hapsira.core import iod
+from hapsira.core.iod import compute_y_vf, tof_equation_y_vf, compute_T_min_gf, find_xy_gf
 from hapsira.iod import izzo
 ```
 
 ## Part 1: Reproducing the original figure
 
-```{code-cell}
+```{code-cell} ipython3
 x = np.linspace(-1, 2, num=1000)
 M_list = 0, 1, 2, 3
 ll_list = 1, 0.9, 0.7, 0, -0.7, -0.9, -1
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 fig, ax = plt.subplots(figsize=(10, 8))
 ax.set_prop_cycle(
     cycler("linestyle", ["-", "--"])
@@ -43,8 +43,8 @@ for M in M_list:
     for ll in ll_list:
         T_x0 = np.zeros_like(x)
         for ii in range(len(x)):
-            y = iod._compute_y(x[ii], ll)
-            T_x0[ii] = iod._tof_equation_y(x[ii], y, 0.0, ll, M)
+            y = compute_y_vf(x[ii], ll)
+            T_x0[ii] = tof_equation_y_vf(x[ii], y, 0.0, ll, M)
         if M == 0 and ll == 1:
             T_x0[x > 0] = np.nan
         elif M > 0:
@@ -86,12 +86,12 @@ ax.set_ylabel("$T$")
 
 ## Part 2: Locating $T_{min}$
 
-```{code-cell}
+```{code-cell} ipython3
 :tags: [nbsphinx-thumbnail]
 
 for M in M_list:
     for ll in ll_list:
-        x_T_min, T_min = iod._compute_T_min(ll, M, 10, 1e-8)
+        x_T_min, T_min = compute_T_min_gf(ll, M, 10, 1e-8)
         ax.plot(x_T_min, T_min, "kx", mew=2)
 
 fig
@@ -99,17 +99,21 @@ fig
 
 ## Part 3: Try out solution
 
-```{code-cell}
+```{code-cell} ipython3
 T_ref = 1
 ll_ref = 0
 
-x_ref, _ = iod._find_xy(
-    ll_ref, T_ref, M=0, numiter=10, lowpath=True, rtol=1e-8
+x_ref, _ = find_xy_gf(
+    ll_ref, T_ref, 
+    0,  # M
+    10, # numiter 
+    True, # lowpath
+    1e-8, # rtol
 )
 x_ref
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 ax.plot(x_ref, T_ref, "o", mew=2, mec="red", mfc="none")
 
 fig
@@ -117,7 +121,7 @@ fig
 
 ## Part 4: Run some examples
 
-```{code-cell}
+```{code-cell} ipython3
 from astropy import units as u
 
 from hapsira.bodies import Earth
@@ -125,7 +129,7 @@ from hapsira.bodies import Earth
 
 ### Single revolution
 
-```{code-cell}
+```{code-cell} ipython3
 k = Earth.k
 r0 = [15945.34, 0.0, 0.0] * u.km
 r = [12214.83399, 10249.46731, 0.0] * u.km
@@ -138,7 +142,7 @@ v0, v = izzo.lambert(k, r0, r, tof)
 v
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 k = Earth.k
 r0 = [5000.0, 10000.0, 2100.0] * u.km
 r = [-14600.0, 2500.0, 7000.0] * u.km
@@ -153,7 +157,7 @@ v
 
 ### Multiple revolutions
 
-```{code-cell}
+```{code-cell} ipython3
 k = Earth.k
 r0 = [22592.145603, -1599.915239, -19783.950506] * u.km
 r = [1922.067697, 4054.157051, -8925.727465] * u.km
@@ -169,20 +173,20 @@ expected_va_r = [-2.45759553, 1.16945801, 0.43161258] * u.km / u.s
 expected_vb_r = [-5.53841370, 0.01822220, 5.49641054] * u.km / u.s
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 v0, v = izzo.lambert(k, r0, r, tof, M=0)
 v
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 _, v_l = izzo.lambert(k, r0, r, tof, M=1, lowpath=True)
 _, v_r = izzo.lambert(k, r0, r, tof, M=1, lowpath=False)
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 v_l
 ```
 
-```{code-cell}
+```{code-cell} ipython3
 v_r
 ```