pythonで微分方程式を解く

前提・実現したいこと

pythonで微分方程式を解きたい
dX/dt=Q*X
Qは６×６の行列で、時間変化しません。

scilabのコードをpythonに書き換えたいです

このサイトを参考にしました
(https://ir.kagoshima-u.ac.jp/?action=repository_action_common_download&item_id=4578&item_no=1&attribute_id=16&file_no=1)

発生している問題・エラーメッセージ

Scilabでは正しいグラフが出るのですがpythonでは、
エラーは出ていませんが明らかに間違っているグラフが表示されます。
[scilabでのグラフ]
[pythonでのグラフ]

該当のソースコード

python
1from re import X
2from control import matlab
3from scipy.integrate.quadpack import quad_explain
4from sympy import *
5import numpy as np
6from control import  ss
7from control.matlab import  initial
8from scipy.linalg.basic import det
9
10import matplotlib.pyplot as plt
11import math
12from control.matlab import acker, initial
13from control import ss
14from scipy.integrate import odeint
15t1 = 100
16th1 = math.pi/4
17th2 = math.pi/6
18x = -0.5
19v = 0
20w1 = 0
21w2 = 0
22g = 9.80665
23M = 1.5#台車の質量
24m1 = 0.15#振り子1の質量
25m2 = 0.15#振り子2の質量
26l1 = 1.2#振り子1の端点と重心の距離
27l2 = 1.2#振り子2の端点と重心の距離
28ux = 0.04#台車の粘性摩擦係数
29uth = 0.05#振り子接続部の粘性摩擦係数
30s = 4*l1*l1*l2*l2*((4*m1+3*m2)*M+(m1+m2)*m1)*m2
31a42 = -6*l1*l1*l2*l2*m2*(2*m1+m2)*(m1+2*m2)*g / s
32a52 = 3*l1*l2*l2*m2*(4*M+4*m1+m2)*(m1+2*m2)*g / s
33a62 = -9*l1*l1*l2*m2*(2*M+m1)*(m1+2*m2)*g / s
34a43 = 6*l1*l1*l2*l2*m1*m2*m2*g / s
35a53 = -9*l1*l2*l2*(2*M+m1)*m2*m2*g / s
36a63 = 3*l1*l1*l2*(4*M*(m1+3*m2)+m1*(m1+4*m2))*m2*g / s
37a44 = -4*l1*l1*l2*l2*(4*m1+3*m2)*m2*ux /s
38a54 = 6*l1*l2*l2*(2*m1+m2)*m2*ux / s
39a64 =-6*l1*l1*l2*m1*m2*ux / s
40a45 = 6*l1*l2*l2*(2*m1+m2)*m2*uth / s
41a55 = -3*l2*l2*(4*M+4*m1+m2)*m2*uth / s
42a65 = 9*l1*l2*(2*M+m1)*m2*uth / s
43a46 = -6*l1*l1*l2*m1*m2*uth / s
44a56 = 9*l1*l2*(2*M+m1)*m2*uth / s
45a66 = -3*l1*l1*(4*M*(m1+3*m2)+m1*(m1+4*m2))*uth / s
46b4 = 4*l1*l1*l2*l2*(4*m1+3*m2)*m2 / s
47b5 = -6*l1*l2*l2*(2*m1+m2)*m2 / s
48b6 = 6*l1*l1*l2*m1*m2 /s
49
50A= np.array([[0,0,0,1,0,0],
51             [0,0,0,0,1,0],
52             [0,0,0,0,0,1],
53             [0,a42,a43,a44,a45,a46],
54             [0,a52,a53,a54,a55,a56],
55             [0,a62,a63,a64,a65,a66]])
56B = np.array([[0],[0],[0],[b4],[b5],[b6]])
57
58C = np.array([[1,0,0,0,0,0],
59              [0,1,0,0,0,0],
60              [0,0,1,0,0,0]])
61
62D = 0
63
64Pole = [-1+1j, -1-1j ,-2+1j,-2-1j,-2,-5] 
65k = matlab.place(A, B, Pole)
66
67Q=A-B*k
68
69def ode(Q,X,t):
70    dXdt=Q*X
71    return dXdt
72
73y0=[x,th1,th2,v,w1,w2]
74
75dt = 0.01
76T = 10.0
77times = np.arange(0.0, T, dt)
78args=(Q,)
79X = odeint(ode, y0, times, args)
80
81plt.plot(times,X[:,0]) 
82plt.grid()
83plt.show()
84plt.figure()

scilab
1clear();
2t1 = 100;
3th1 = %pi/4;
4th2 = -%pi/6;
5x = -0.5;
6v = 0;
7w1 = 0;
8w2 = 0;
9g = 9.80665;
10M = 1.5;//台車の質量
11m1 = 0.15;//振り子1の質量
12m2 = 0.15;//振り子2の質量
13l1 = 1.2;//振り子1の端点と重心の距離
14l2 = 1.2;//振り子2の端点と重心の距離
15ux = 0.04;//台車の粘性摩擦係数
16uth = 0.05;//振り子接続部の粘性摩擦係数
17s = 4*l1*l1*l2*l2*((4*m1+3*m2)*M+(m1+m2)*m1)*m2;
18a42 = -6*l1*l1*l2*l2*m2*(2*m1+m2)*(m1+2*m2)*g / s;
19a52 = 3*l1*l2*l2*m2*(4*M+4*m1+m2)*(m1+2*m2)*g / s;
20a62 = -9*l1*l1*l2*m2*(2*M+m1)*(m1+2*m2)*g / s;
21a43 = 6*l1*l1*l2*l2*m1*m2*m2*g / s;
22a53 = -9*l1*l2*l2*(2*M+m1)*m2*m2*g / s;
23a63 = 3*l1*l1*l2*(4*M*(m1+3*m2)+m1*(m1+4*m2))*m2*g / s;
24a44 = -4*l1*l1*l2*l2*(4*m1+3*m2)*m2*ux / s;
25a54 = 6*l1*l2*l2*(2*m1+m2)*m2*ux / s;
26a64 =-6*l1*l1*l2*m1*m2*ux / s;
27a45 = 6*l1*l2*l2*(2*m1+m2)*m2*uth / s;
28a55 = -3*l2*l2*(4*M+4*m1+m2)*m2*uth / s;
29a65 = 9*l1*l2*(2*M+m1)*m2*uth / s;
30a46 = -6*l1*l1*l2*m1*m2*uth / s;
31a56 = 9*l1*l2*(2*M+m1)*m2*uth / s;
32a66 = -3*l1*l1*(4*M*(m1+3*m2)+m1*(m1+4*m2))*uth / s;
33b4 = 4*l1*l1*l2*l2*(4*m1+3*m2)*m2 / s;
34b5 = -6*l1*l2*l2*(2*m1+m2)*m2 / s;
35b6 = 6*l1*l1*l2*m1*m2 /s;
36A = [0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;0 a42 a43 a44 a45 a46;0 a52 a53 a54 a55 a56;0 a62 a63 a64 a65 a66];
37B = [0 ; 0 ; 0 ; b4 ; b5 ; b6];
38C = [1 0 0 0 0 0;0 1 0 0 0 0 ;0 0 1 0 0 0];
39D = 0;
40poles = [-1+1*%i,-1-1*%i,-2+1*%i,-2-1*%i,-2,-5];//極
41f = ppol(A,B,poles);
42function Xdot=F(t,X)
43    Xdot=(A-B*f)*X
44endfunction
45X0=[x;th1;th2;v;w1;w2];//初期値
46X=X0;
47t0=0;
48t=0:0.1:10;
49X=ode(X0,t0,t,F);//
50
51xgrid();
52plot(t,X(1,:))
53xlabel("time[s]",'FontSize',4);
54ylabel("x[m]",'FontSize',4);
55h1=legend(['x']);
56

試したこと

def の中身が間違っていると思うのですが、どのように書けばよいでしょうか？
どうかご教授お願い致します。

python
1Q=A-B*k
2
3def ode(Q,X,t):
4    dXdt=Q*X
5    return dXdt
6
7y0=[x,th1,th2,v,w1,w2]
8
9dt = 0.01
10T = 10.0
11times = np.arange(0.0, T, dt)
12args=(Q,)
13X = odeint(ode, y0, times, args)
14
15plt.plot(times,X[:,0])
16