Solving second order ODE numerically

Here is how I solved it in MATLAB:

h0 = 1;  % initial value for h
dhdx0 = 1e-5;  % initial value for dh/dx
 
xspan = [0 1e4];
y0 = [h0 dhdx0];
 
c = (5/3) * (1e-3)^2;
 
[x,y] = ode45(@(x,y) odefcn(x,y,c), xspan, y0);
 
plot(x,y(:,1),'-o',x,y(:,2),'-.')
 
function dydt = odefcn(x,y,c)
    dydt = zeros(2,1);
    dydt(1) = y(2);
    dydt(2) = c * (y(1)^(-3) - y(1)^(-6));
end

I tried to do the same in Python, but get strange results (think that MATLABs smart mesh is the difference):

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

h0 = 1  # initial value for h
dhdx0 = 1e-5  # initial value for dh/dx

m = GEKKO(remote=False)
m.time = np.linspace(0, 1e4, 250)
y = m.Var(h0)
dy = m.Var(dhdx0)

c = (5/3) * (1e-3)**2
m.Equation(y.dt() == dy)
m.Equation(dy.dt() == c * (y**(-3) - y**(-6)))

m.options.IMODE = 6
m.options.NODES = 3
m.solve(disp=False)

plt.plot(m.time, y.value, '.-')
plt.show()

EDIT: as suggested by @LutzLehmann, the equivalent to ODE45 in Python is scipy.integrate.solve_ivp:

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

h0 = 1  # initial value for h
dhdx0 = 1e-5  # initial value for dh/dx


def f(x, y, c):
    return y[1], c * (y[0]**(-3) - y[0]**(-6))


xspan = np.linspace(0, 1e4, 100)
c = (5/3) * (1e-3)**2
sol = solve_ivp(f, [xspan[0], xspan[-1]], [h0, dhdx0],
                args=(c,), t_eval=xspan, rtol=1e-5)

plt.plot(sol.t, sol.y[0], '.-', label='Output')
plt.show()