How to plot figures side by side in custom function?

Question

I have a few plots using Matplotlib and partial_dependence that I use the below code to generate:

def plot_pdp(model, X, feature, target=False, return_pd=False, y_pct=True, figsize=(10,9), norm_hist=True, dec=.5):
    # Get partial dependence
    pardep = partial_dependence(model, X, [feature])
    
    # Get min & max values
    xmin = pardep[1][0].min()
    xmax = pardep[1][0].max()
    ymin = pardep[0][0].min()
    ymax = pardep[0][0].max()
    
    # Create figure
    fig, ax1 = plt.subplots(figsize=figsize)
    ax1.grid(alpha=.5, linewidth=1)
    
    # Plot partial dependence
    color = 'tab:blue'
    ax1.plot(pardep[1][0], pardep[0][0], color=color)
    ax1.tick_params(axis='y', labelcolor=color)
    ax1.set_xlabel(feature, fontsize=14)
    
    tar_ylabel = ': {}'.format(target) if target else ''
    ax1.set_ylabel('Partial Dependence{}'.format(tar_ylabel), color=color, fontsize=14)
    
    tar_title = target if target else 'Target Variable'
    ax1.set_title('Relationship Between {} and {}'.format(feature, tar_title), fontsize=16)
    
    if y_pct and ymin>=0 and ymax<=1:
        # Display yticks on ax1 as percentages
        fig.canvas.draw()
        labels = [item.get_text() for item in ax1.get_yticklabels()]
        labels = [int(np.float(label)*100) for label in labels]
        labels = ['{}%'.format(label) for label in labels]
        ax1.set_yticklabels(labels)
    
    # Plot line for decision boundary
    #ax1.hlines(dec, xmin=xmin, xmax=xmax, color='black', linewidth=2, linestyle='--', label='Decision Boundary')
    #ax1.legend()

    ax2 = ax1.twinx()
    color = 'tab:red'
    ax2.hist(X[feature], bins=80, range=(xmin, xmax), alpha=.25, color=color, density=norm_hist)
    ax2.tick_params(axis='y', labelcolor=color)
    ax2.set_ylabel('Distribution', color=color, fontsize=14)
    
    if y_pct and norm_hist:
        # Display yticks on ax2 as percentages
        fig.canvas.draw()
        labels = [item.get_text() for item in ax2.get_yticklabels()]
        labels = [int(np.float(label)*100) for label in labels]
        labels = ['{}%'.format(label) for label in labels]
        ax2.set_yticklabels(labels)

    plt.show()
    
    if return_pd:
        return pardep

Then I can get a plot for each column using :

for i in X.columns:
    plot_pdp(model, X, i, target='nof_visit')

In a Jupyter notebook the plots are shown one below another. How can I change it so that it is shown side by side?

Edit: When I try subplots on this line 'fig, ax1 = plt.subplots(1, len(X.columns) ,figsize=figsize)', I get the below error:

AttributeError                            Traceback (most recent call last)
Input In [174], in <cell line: 71>()
     69         return pardep
     71 for i in X.columns:
---> 72     plot_pdp(model, X, i, target='nof_visit')

Input In [174], in plot_pdp(model, X, feature, target, return_pd, y_pct, figsize, norm_hist, dec)
     11     # Create figure
     12     fig, ax1 = plt.subplots(1, len(X.columns) ,figsize=figsize)
---> 13     ax1.grid(alpha=.5, linewidth=1)
     15 #     SMALL_SIZE = 20
     16 #     MEDIUM_SIZE = 40
     17 #     BIGGER_SIZE = 40
   (...)
     27     
     28     # Plot partial dependence
     29     color = 'tab:blue'

AttributeError: 'numpy.ndarray' object has no attribute 'grid'