How do I get n_samples in one call of pm.sample_posterior_predictive() in pymc?

I want to generate multiple samples from one function call of pymc's sample_posterior_predictive(). In the previous version pymc3 there was an argument called 'samples', here is an example https://www.pymc.io/projects/examples/en/latest/gaussian_processes/GP-Marginal.html. However this argument is not part of the API anymore in pymc version '5.1.2'.

I am further using matplotlib version '3.7.1', numpy version '1.24.2' and python version '3.11.0'.

Here is a minimal example (as minimal as possible) showing a workaround, however I really want to avoid the for-loop and generate n_samples from one function call.

import matplotlib.pyplot as plt
import numpy as np
import pymc as pm

# set the seed
np.random.seed(1)
# create random data
n = 50  # The number of data points
X = np.linspace(0, np.pi*2, n)[:, None]  # The inputs to the GP, they must be arranged as a column vector
y = 2*np.sin(0.25*2*np.pi*X[:,0])*X[:,0] + 2

# setup model
with pm.Model() as model:
    ℓ = pm.Gamma("ℓ", alpha=2, beta=1)
    η = pm.HalfCauchy("η", beta=5)

    cov = η**2 * pm.gp.cov.Matern52(1, ℓ)
    gp = pm.gp.Marginal(cov_func=cov)

    σ = pm.HalfCauchy("σ", beta=5)
    y_ = gp.marginal_likelihood("y", X=X, y=y, sigma=σ)

    mp = pm.find_MAP()
    
# new values from x=0 to x=20
X_new = np.linspace(0, 20, 600)[:, None]

# add the GP conditional to the model, given the new X values
with model:
    f_pred = gp.conditional("f_pred", X_new)
    
########################################################################################
# How do I get (n_samples, 600) in one call of pm.sample_posterior_predictive()?
n_samples = 4
sample_list = []
for i in range(n_samples):
    with model:
        pred_samples = pm.sample_posterior_predictive([mp], var_names=["f_pred"])
        sample_list.append(pred_samples)
########################################################################################

# plot result
fig, ax = plt.subplots( figsize=(12, 5))
for pred_samples in sample_list:
    f_pred = pred_samples.posterior_predictive["f_pred"].sel(chain=0)
    ax.plot(X_new[:,0], f_pred[0,:],  alpha=0.1, color = 'blue')
    
# plot the data and the true latent function
ax.plot(X, y, "ok", ms=3, alpha=0.5, label="Observed data")
plt.show()