diff --git a/floris/core/core.py b/floris/core/core.py
index 857b90fa4..afe4d9c76 100644
--- a/floris/core/core.py
+++ b/floris/core/core.py
@@ -264,8 +264,6 @@ def solve_for_points(self, x, y, z):
         else:
             full_flow_sequential_solver(self.farm, self.flow_field, field_grid, self.wake)
 
-        return self.flow_field.u_sorted[:,:,0,0] # Remove turbine grid dimensions
-
     def solve_for_velocity_deficit_profiles(
         self,
         direction: str,
@@ -319,7 +317,8 @@ def solve_for_velocity_deficit_profiles(
         y = np.squeeze(y, axis=0) + y_start
         z = np.squeeze(z, axis=0) + reference_height
 
-        u = self.solve_for_points(x.flatten(), y.flatten(), z.flatten())
+        self.solve_for_points(x.flatten(), y.flatten(), z.flatten())
+        u = self.flow_field.u_sorted[:, :, 0, 0]
         u = np.reshape(u[0, :], (n_lines, resolution))
         velocity_deficit = (homogeneous_wind_speed - u) / homogeneous_wind_speed
 
diff --git a/floris/floris_model.py b/floris/floris_model.py
index aaa384029..d944b6209 100644
--- a/floris/floris_model.py
+++ b/floris/floris_model.py
@@ -1368,7 +1368,11 @@ def get_plane_of_points(
 
         return df
 
-    def sample_flow_at_points(self, x: NDArrayFloat, y: NDArrayFloat, z: NDArrayFloat):
+    def sample_flow_at_points(self,
+                              x: NDArrayFloat,
+                              y: NDArrayFloat,
+                              z: NDArrayFloat,
+                              run_solver: bool = True):
         """
         Extract the wind speed at points in the flow.
 
@@ -1376,6 +1380,7 @@ def sample_flow_at_points(self, x: NDArrayFloat, y: NDArrayFloat, z: NDArrayFloa
             x (1DArrayFloat | list): x-locations of points where flow is desired.
             y (1DArrayFloat | list): y-locations of points where flow is desired.
             z (1DArrayFloat | list): z-locations of points where flow is desired.
+            run_solver (bool): If `True` (default), run the solver.
 
         Returns:
             3DArrayFloat containing wind speed with dimensions
@@ -1386,9 +1391,17 @@ def sample_flow_at_points(self, x: NDArrayFloat, y: NDArrayFloat, z: NDArrayFloa
         if not len(x) == len(y) == len(z):
             raise ValueError("x, y, and z must be the same size")
 
-        return self.core.solve_for_points(x, y, z)
+        if run_solver:
+            self.core.solve_for_points(x, y, z)
 
-    def sample_ti_at_points(self, x: NDArrayFloat, y: NDArrayFloat, z: NDArrayFloat):
+        # Remove grid dimensions and return sorted wind speed field.
+        return self.core.flow_field.u_sorted[:, :, 0, 0]
+
+    def sample_ti_at_points(self,
+                            x: NDArrayFloat,
+                            y: NDArrayFloat,
+                            z: NDArrayFloat,
+                            run_solver: bool = True):
         """
         Extract the turbulence intensity at points in the flow.
 
@@ -1396,13 +1409,14 @@ def sample_ti_at_points(self, x: NDArrayFloat, y: NDArrayFloat, z: NDArrayFloat)
             x (1DArrayFloat | list): x-locations of points where TI is desired.
             y (1DArrayFloat | list): y-locations of points where TI is desired.
             z (1DArrayFloat | list): z-locations of points where TI is desired.
+            run_solver (bool): If `True` (default), run the solver.
 
         Returns:
             3DArrayFloat containing turbulence intensity with dimensions
             (# of findex, # of sample points)
         """
-
-        self.sample_flow_at_points(x, y, z) # Solve, but ignore returned velocities
+        if run_solver:
+            self.core.solve_for_points(x, y, z)
 
         # Remove grid dimensions and return sorted TI field
         return self.core.flow_field.turbulence_intensity_field_sorted[:, :, 0, 0]