Issue with GEKKO bspline for coolprop data in two-phase region

I am trying to use coolprop inside my dynamic thermal simulation in gekko. I learned that rather than implementing coolprop calls directly bspline or cspline objects have to be used to allow gekko to find the derivatives. However, when attempting to do so in the two-phase region of the fluid, the bspline object does not lead to the correct results. Below is a minimal working example where I show the gekko results compared to a 2d spline from scipy. I am using a spline calculating the fluid´s temperature depending on the pressure and specific enthalpy.

import numpy as np
from CoolProp.CoolProp import PropsSI as PS
from gekko import GEKKO
from matplotlib import pyplot as plt
from scipy.interpolate import bisplrep, bisplev


def coolprop_bspline():
    # Gekko Setup
    m = GEKKO(remote=False)
    t_steps = 90
    m.time = np.linspace(0, 90, t_steps)

    # Initial values
    xi = 0.0
    p_i = 10e5
    tau_i = PS('T', 'Q', xi, 'P', p_i, 'R134a')
    h_i = PS('H', 'Q', xi, 'P', p_i, 'R134a')
    tau_low = tau_i - 30
    tau_high = tau_i + 30

    # Prepare spline data
    h_low = PS('H', 'T', tau_low, 'P', p_i, 'R134a')
    h_high = PS('H', 'T', tau_high, 'P', p_i, 'R134a')
    h_range = np.squeeze(np.linspace(h_low, h_high, 50))
    p_range = np.linspace(9.0e5, 11e5, 10)
    h_mesh, p_mesh = np.meshgrid(h_range, p_range, indexing='ij')
    tau_mesh = np.reshape(
        PS('T', 'H', np.reshape(h_mesh, (h_mesh.size,)), 'P', np.reshape(p_mesh, (p_mesh.size,)), 'R134a'),
        h_mesh.shape)

    # Variables
    h = m.Var(h_i)
    p = m.Var(p_i)
    tau = m.Var(tau_i)

    # Equations
    m.Equations(
        [
            h.dt() == 1000,
            p.dt() == 0,
        ]
    )
    m.bspline(h, p, tau, h_range, p_range, tau_mesh, data=True, kx=3, ky=3)

    # Execution
    m.options.IMODE = 4
    m.options.SOLVER = 3
    m.solve(disp=True)

    # Evaluate comparative spline from Scipy
    spl = bisplrep(h_mesh, p_mesh, tau_mesh, kx=3, ky=3)
    h_new = np.linspace(h_low, h_high, 200)
    p_new = np.linspace(9.0e5, 11e5, 200)
    h_new, p_new = np.meshgrid(h_new, p_new, indexing='ij')
    tau_interp = bisplev(h_new[:, 0], p_new[0, :], spl)

    # Plotting

    # Plot result from Gekko
    plt.figure(1)
    plt.plot(h.value,tau.value)
    plt.grid()
    plt.xlabel('Enthalpy [J/kg]')
    plt.ylabel('Temperature [K]')
    plt.title('Gekko results with p = 1e6')

    # Plot comparative spline from scipy
    fig = plt.figure(2)
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(h_new, p_new, tau_interp)
    ax.view_init(20, -100)
    plt.xlabel('Enthalpy [J/kg]')
    plt.ylabel('Pressure [Pa]')
    ax.set_zlabel('Temperature [K]')
    plt.title('Correct BSpline')

    # Plot original data from coolprop
    fig = plt.figure(3)
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(h_mesh, p_mesh, tau_mesh)
    ax.view_init(20,-100)
    plt.xlabel('Enthalpy [J/kg]')
    plt.ylabel('Pressure [Pa]')
    ax.set_zlabel('Temperature [K]')
    plt.title('Raw coolprop data')

    plt.show()


if __name__ == '__main__':
    coolprop_bspline()

Below are the resulting graphs. When comparing the Gekko result with the other two, it is obvious that the temperature should not suddenly increase to the maximum value covered by the splines at an enthalpy of around 335000 J/kg.

Does anyone have an idea of how to adjust the bspline object to correctly fit the data? So far I have already tried different polynomial degrees (kx, ky) as well as logarithmic inputs for pressure and enthalpy.