diff --git a/lumicks/pylake/fitting/model.py b/lumicks/pylake/fitting/model.py
index 6bff76278..1701c6ac1 100644
--- a/lumicks/pylake/fitting/model.py
+++ b/lumicks/pylake/fitting/model.py
@@ -1,5 +1,6 @@
 import uuid
 import inspect
+import warnings
 from copy import deepcopy
 from collections import OrderedDict
 
@@ -241,7 +242,12 @@ def __repr__(self):
         return model_info
 
     def invert(
-        self, independent_min=0.0, independent_max=np.inf, interpolate=False, independent_step=1e-2
+        self,
+        independent_min=None,
+        independent_max=None,
+        interpolate=False,
+        method="exact",
+        independent_step=1e-2,
     ):
         """Invert this model.
 
@@ -264,10 +270,46 @@ def invert(
             `interpolate` is set to `True`.
         interpolate : bool, optional
             Use interpolation method rather than numerical inversion.
+        method : {"exact", "interp_root", "interp_lsq"}
+            - "exact" : Invert each point separately (exact, but slow).
+            - "interp_root" : Interpolate forward model as lookup table for inversion. Appropriate range is selected by rootfinding.
+            - "interp_lsq" : Interpolate forward model as lookup table for inversion. Appropriate range is selected by least_squares (deprecated).
+              Can only be used when track is equidistantly sampled.
+
         independent_step : float, optional
             Interpolation step size. Default: 1e-2
         """
-        return InverseModel(self, independent_min, independent_max, interpolate, independent_step)
+        if method not in ("exact", "interp_root", "interp_lsq"):
+            raise RuntimeError("Interpolation method should either be rootfinding or least_squares")
+
+        if interpolate:
+            method = "exact"
+            warnings.warn(
+                RuntimeWarning("The parameter `interpolate` is deprecated. Use method parameter.")
+            )
+
+        if method == "exact":
+            return InverseModel(
+                self, independent_min, independent_max, False, independent_step, rootfinding=False
+            )
+        elif method == "interp_root":
+            return InverseModel(
+                self,
+                1e-4 if independent_min is None else independent_min,
+                1e5 if independent_max is None else independent_max,
+                True,
+                independent_step,
+                rootfinding=True,
+            )
+        elif method == "interp_lsq":
+            return InverseModel(
+                self,
+                0.0 if independent_min is None else independent_min,
+                np.inf if independent_max is None else independent_max,
+                True,
+                independent_step,
+                rootfinding=False,
+            )
 
     def subtract_independent_offset(self):
         """
@@ -611,6 +653,7 @@ def __init__(
         independent_max=np.inf,
         interpolate=False,
         independent_step=1e-2,
+        rootfinding=False,
     ):
         """
         Combine two model outputs to form a new model (addition).
@@ -627,6 +670,8 @@ def __init__(
             Use interpolation approximation. Default: False.
         independent_step : float, optional
             Interpolation step size. Default: 1e-2
+        rootfinding : bool
+            Use rootfinding rather than least squares to invert model.
 
         Raises
         ------
@@ -641,6 +686,7 @@ def __init__(
         self.independent_min = independent_min
         self.independent_max = independent_max
         self.independent_step = independent_step
+        self.rootfinding = rootfinding
 
     @property
     def dependent(self):
@@ -678,6 +724,7 @@ def _raw_call(self, independent, param_vector):
                 lambda f_trial: self.model._raw_call(f_trial, param_vector),
                 lambda f_trial: self.model.derivative(f_trial, param_vector),
                 dx=self.independent_step,
+                rootfinding=self.rootfinding,
             )
         else:
             return invert_function(
@@ -687,6 +734,7 @@ def _raw_call(self, independent, param_vector):
                 self.independent_max,
                 lambda f_trial: self.model._raw_call(f_trial, param_vector),  # Forward model
                 lambda f_trial: self.model.derivative(f_trial, param_vector),
+                rootfinding=self.rootfinding,
             )
 
     @property