diff --git a/docs/examples/shading/plot_partial_module_shading_simple.py b/docs/examples/shading/plot_partial_module_shading_simple.py
index ce031eff10..4d04926e8d 100644
--- a/docs/examples/shading/plot_partial_module_shading_simple.py
+++ b/docs/examples/shading/plot_partial_module_shading_simple.py
@@ -38,7 +38,7 @@
 from pvlib import pvsystem, singlediode
 import pandas as pd
 import numpy as np
-from scipy.interpolate import interp1d
+from scipy.interpolate import make_interp_spline
 import matplotlib.pyplot as plt
 
 from scipy.constants import e as qe, k as kB
@@ -178,10 +178,13 @@ def plot_curves(dfs, labels, title):
 
 
 def interpolate(df, i):
-    """convenience wrapper around scipy.interpolate.interp1d"""
-    f_interp = interp1d(np.flipud(df['i']), np.flipud(df['v']), kind='linear',
-                        fill_value='extrapolate')
-    return f_interp(i)
+    x = np.flipud(df['i'])
+    y = np.flipud(df['v'])
+    
+    # Create a spline interpolation with linear (k=1) and extrapolation (default behavior)
+    spline = make_interp_spline(x, y, k=1, bc_type='clamped')  # Extrapolation is handled by default
+    
+    return spline(i)
 
 
 def combine_series(dfs):
diff --git a/pvlib/iam.py b/pvlib/iam.py
index 161de84589..ba03dd60b7 100644
--- a/pvlib/iam.py
+++ b/pvlib/iam.py
@@ -470,7 +470,8 @@ def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
     '''
     # Contributed by Anton Driesse (@adriesse), PV Performance Labs. July, 2019
 
-    from scipy.interpolate import interp1d
+    from scipy.interpolate import make_interp_spline
+    import numpy as np
 
     # Scipy doesn't give the clearest feedback, so check number of points here.
     MIN_REF_VALS = {'linear': 2, 'quadratic': 3, 'cubic': 4, 1: 2, 2: 3, 3: 4}
@@ -483,8 +484,10 @@ def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
         raise ValueError("Negative value(s) found in 'iam_ref'. "
                          "This is not physically possible.")
 
-    interpolator = interp1d(theta_ref, iam_ref, kind=method,
-                            fill_value='extrapolate')
+    theta_ref = np.asarray(theta_ref)
+    iam_ref = np.asarray(iam_ref)
+
+    interpolator = make_interp_spline(theta_ref, iam_ref, k=method, bc_type='clamped')
     aoi_input = aoi
 
     aoi = np.asanyarray(aoi)
diff --git a/pvlib/spectrum/response.py b/pvlib/spectrum/response.py
index 4da92bb32a..2e374c1882 100644
--- a/pvlib/spectrum/response.py
+++ b/pvlib/spectrum/response.py
@@ -6,7 +6,7 @@
 import numpy as np
 import pandas as pd
 import scipy.constants
-from scipy.interpolate import interp1d
+from scipy.interpolate import make_interp_spline
 
 
 _PLANCK_BY_LIGHT_SPEED_OVER_ELEMENTAL_CHARGE_BY_BILLION = (
@@ -67,12 +67,11 @@ def get_example_spectral_response(wavelength=None):
         resolution = 5.0
         wavelength = np.arange(280, 1200 + resolution, resolution)
 
-    interpolator = interp1d(SR_DATA[0], SR_DATA[1],
-                            kind='cubic',
-                            bounds_error=False,
-                            fill_value=0.0,
-                            copy=False,
-                            assume_sorted=True)
+    x_data = np.sort(SR_DATA[0])
+    y_data = SR_DATA[1]
+
+    interpolator = make_interp_spline(x_data, y_data, k=3, bc_type='clamped')
+
 
     sr = pd.Series(data=interpolator(wavelength), index=wavelength)