#### Differences between Matlab ode45 and Scipy odeint: same model different results

EDIT

So I was able to produce the correct plot, but it was only after adjusting the sampling interval to the following:

`t_list = np.linspace(0, 30, 100)`

which prints X as:

``````[1.         1.         1.         1.         1.         1.
1.         1.         1.         1.         1.         1.
1.         1.         1.         1.         1.         1.
0.91265299 0.8107059  0.7366542  0.68370578 0.64633005 0.62021062
0.6020953  0.58960093 0.58101775 ...]
``````

But that begs the question: why is this system so dependent on the sampling intervals?

END EDIT

I am trying to recreate a simple matlab system of differential equation into python using scipy. I do not know why I am receiving a `RuntimeWarning: invalid value encountered in double_scalars` in the scipy execution. Am I missing an optional parameter in the odeint call?

python:

``````import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def model(y0, t):
x = y0
y = y0
z = y0
if t <= 10:
sys_input = 1.0
else:
sys_input = 0.75
a = 1.0
b = 1.0
c = 1.0
E = 1.0
dxdt = sys_input - a * E * (x ** 0.5)
dydt = a * E * (x ** 0.5) - b * (y ** 0.5)
dzdt = b * (y ** 0.5) - c * (z ** 0.5)
return [dxdt, dydt, dzdt]

t_list = np.linspace(0, 30, 31)

# Initial conditions vector
yi = [1.0, 1.0, 1.0]
ret = odeint(model, y0=yi, t=t_list)
X = ret[:, 0]
print(X)
``````

and it prints the following:

``````<input>:18: RuntimeWarning: invalid value encountered in double_scalars
<input>:19: RuntimeWarning: invalid value encountered in double_scalars
[ 1.  1.  1.  1. nan nan nan nan nan nan nan nan nan nan nan nan nan nan
nan nan nan nan nan nan nan nan nan nan nan nan nan]
``````

where as the following matlab code produces a continuous result:

``````tspan = [0,30];
x0 = 1.0; % Initial value of x
y0 = 1.0; % Initial value of y
z0 = 1.0; % Initial value of z
initial_values = [x0; y0; z0]; % Initial value of the vector w
[T,R] = ode45(@(t,w) diff_eq(t,w),tspan,initial_values);
X = R(:,1);
Y = R(:,2);
Z = R(:,3);

for i = 1: length(X)
if(mod(i, 10)==0 &&  i > 1)
disp(' ');
end
fprintf('X[%i] = %.2f, ', i, X(i));
end
disp(' ');

function dw_vectordt = diff_eq(t,w_vector)
x = w_vector(1);
y = w_vector(2);
z = w_vector(3);
if (t<=10)
sys_input= 1.0;
else
sys_input=0.75;
end
a = 1.0;
b = 1.0;
c = 1.0;
E = 1.0;
dxdt = sys_input-a*E*x^(0.5);
dydt = a*E*x^(0.5)-b*y^(0.5);
dzdt = b*y^(0.5)-c*z^(0.5);
dw_vectordt = [dxdt; dydt; dzdt];
end
``````

the print statement returns:

``````X = 1.00, X = 1.00, X = 1.00, X = 1.00, X = 1.00, X = 1.00, X = 1.00, X = 1.00, X = 1.00,
X = 1.00, X = 1.00, X = 1.00, X = 1.00, X = 1.01, X = 0.99, X = 0.92, X = 0.82, X = 0.78, X = 0.74,
X = 0.71, X = 0.68, X = 0.66, X = 0.64, X = 0.63, X = 0.61, X = 0.60, X = 0.59, X = 0.59, X = 0.58,
X = 0.58, X = 0.57, X = 0.57, X = 0.57, X = 0.57, X = 0.57, X = 0.56, X = 0.56, X = 0.56, X = 0.56,
X = 0.56, X = 0.56, X = 0.56, X = 0.56, X = 0.56, X = 0.56, X = 0.56, X = 0.56, X = 0.56, X = 0.56,
X = 0.56, X = 0.56, X = 0.56, X = 0.56,
``````

Source: Python Questions