scipy.integrate.solve_ivp diverges on a state space simulation well finished in MATLAB

Question

I tried to simulate a state space model with MATLAB ode45, then I tried the same work in Python with scipy.integrate.solve_ivp. As it is obviously shown in this post pictures. Python simulation diverges for no good reason. The solvers message is "Required step size is less than spacing between numbers." but adding time steps is not a solution.

Here is the MATLAB code for the time interval of half a second following a link for the plot of 173rd state: [1]: https://i.stack.imgur.com/mMdNQ.png

C_static=csvread('C_static.csv');
M_static=csvread('M_static.csv');
B_static=csvread('B_static.csv');
CY_static=csvread('CY_static.csv');
DY_static=csvread('DY_static.csv');
dynamoterm = csvread('dynamoterm.csv');

C_static=0*C_static;

n2panto=dynamoterm(6,1);
n2cw=dynamoterm(6,2);

k_dynamic = KCdyna(0,dynamoterm);
K_total = K_static;
K_total(n2cw+1:n2panto+1,n2cw+1:n2panto+1)=K_total(n2cw+1:n2panto+1,n2cw+1:n2panto+1)+k_dynamic;
A_static = [0*K_static, eye(length(B_static));
           -M_static\K_static, -M_static\C_static];

Bu = [0*B_static;
      M_static\B_static];
inc0 = -A_static\Bu;

M_in=inv(M_static);
M_cwp=M_in(:,n2cw+1:n2panto+1);

timer=tic;
[T, Y] = ode45(@(t,X) Asol(X,t,A_static,M_cwp,Bu,n2cw,n2panto,dynamoterm),[0,0.5],inc0);
output=[CY_static,0*CY_static]*Y'+DY_static*ones(1,length(T));
figure
plot(T,output(173,:));
stopwatch=toc(timer);

function dx=Asol(X,t,A_static,M_cwp,Bu,n2cw,n2panto,dynamoterm)

[k_dynamic]=KCdyna(t,dynamoterm);
A=A_static;
A(n2panto+4:2*(n2panto+3),n2cw+1:n2panto+1)=A_static(n2panto+4:2*(n2panto+3),n2cw+1:n2panto+1)-M_cwp*k_dynamic;
dx=A*X+Bu;

end

[![MATLAB simulation plot of the 173rd state][1]][1]

Here is my similar work in Python for the time interval of half a second following a link for the plot of 173rd state: [2]: https://i.stack.imgur.com/LOg2j.png

from KCdyna2 import K_dyn
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# Imports matrices via .csv file
M = np.genfromtxt('Excel\dyn\M_static.csv', delimiter=',')
C = np.genfromtxt('Excel\dyn\C_static.csv', delimiter=',')
C = np.zeros(np.shape(C))
K_static = np.genfromtxt('Excel\dyn\K_static.csv', delimiter=',')
B = np.genfromtxt('Excel\dyn\B_static.csv', delimiter=',')
dyn_trm = np.genfromtxt('Excel\dyn\dynamoterm.csv', delimiter=',')


# Slice addresses
n2cw = int(dyn_trm[5, 1]) # Slice beginning
n2panto = int(dyn_trm[5, 0]) # Slice finishing

# Time interval for solution
time_interval = [0, 0.5]
times = np.linspace(time_interval[0], time_interval[1], 50000)

M_inv = np.linalg.inv(M)

K_total = K_static

K_total[n2cw:n2panto + 1,
        n2cw:n2panto + 1] += K_dyn(0, dyn_trm)

# System dynamics matrix
A_static = np.block([[np.zeros((len(M_inv), len(M_inv)), dtype='uint8'), np.eye(len(M_inv), dtype='uint8')],
                 [np.matmul(-M_inv, K_total), np.matmul(-M_inv, C)]])

Bu = np.block([[np.zeros((len(B), 1), dtype='uint8')],
               [np.matmul(M_inv, B).reshape(len(B), 1)]])

inc0 = np.matmul(-np.linalg.inv(A_static), Bu)

M_inv = np.linalg.inv(M)
M_cwp = M_inv[:, n2cw:n2panto + 1]

def SttSpcEq(t, x, M_cwp, A_st, Bu, dynamoterm, n2cw,n2panto):
    K_dynamic = K_dyn(t, dynamoterm)
    A = A_st
    A[n2panto + 3:2*(n2panto+3),
      n2cw:n2panto + 1] -= np.matmul(M_cwp, K_dynamic)
    return (np.matmul(A, x.reshape(len(Bu), 1)) + Bu).reshape(len(Bu), )

soln = solve_ivp(SttSpcEq,
                time_interval,
                inc0.reshape(len(inc0),),
                method='RK45', t_eval=times, args=(M_cwp, A_static, Bu, dyn_trm, n2cw, n2panto))
print(soln.message)
plt.plot(soln.t, soln.y[172])
plt.show() ```

[![Python simulation plot of the 173rd state][2]][2]