scipy: Bug in Solve_ivp with BDF and Radau solvers and numpy arrays

It seems Scipy integration module solve_ivp when used with BDF and Radau solvers appends rounded numbers to the array. This will then cause an error of inconsistent dimensional further on in my code. I wrote a simple example of a system with (2,1)-dimensional state (i.e. a second order system), which prints the state y as it runs in the function f(t,y), at one point, it appends the rounded number and gave me

[[ 13.67999979  13.67999999]
 [-15.99       -15.98999976]]

The second column is unnecessary and these added numbers will cause error in my original code.

Reproducing code example:

from numpy import array, matmul
from scipy.integrate import solve_ivp
from numpy import size
def sys(t,x,A,xeq,matmul,size):
    if size(x)==4:
        print(x)
    xdot=matmul(A,x-xeq)
    return xdot

T=10
A=array([[-1,2],[-2,-6]])
x0=array([13.67999999,-15.989999999])
xeq=array([[0.99999999],[2.999999999]])
fun=lambda t, x: sys(t,x,A,xeq,matmul,size)
sol=solve_ivp(fun,[0,T],x0,method='BDF',vectorized=True)

The output is:

[[ 13.67999979  13.67999999]
 [-15.99       -15.98999976]]

Scipy/Numpy/Python version information:

1.1.0 1.14.4 sys.version_info(major=3, minor=6, micro=5, releaselevel='final', serial=0)