x = np.linspace(-5, 5, 100)

with pm.Model() as model:
    # Set mutable data
    X = pm.MutableData("X", x[:, None])

    m = [5]
    c = 5.0

    X_mu = np.mean(x, axis=0)
    Xs = X - X_mu
    
    ### WITH EVAL, NO BUG.  WITHOUT EVAL, BUG
    L = pm.gp.hsgp_approx.set_boundary(Xs, c).eval()
    
    eigvals = pm.gp.hsgp_approx.calc_eigenvalues(L, m, tl=pt)
    phi = pm.gp.hsgp_approx.calc_eigenvectors(Xs, L, eigvals, m, tl=pt)

phi1 = phi.eval()

x_new = np.linspace(-5, 10, 200)

with model:
    pm.set_data({"X": x_new[:, None]})
    
plt.plot(x, phi1, 'k');
plt.plot(x_new, phi.eval(), 'b', alpha=0.5);