Solving two coupled non-linear second order differential equations numerically using python

Please help with the following Python code for a nonlinear coupled oscillator, I am getting an error, I don’t understand the problem.

from scipy.optimize import minimize
from scipy.integrate import odeint
from qutip import *
import scipy.special as sp
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

def vectorfield(w, t, p):
    """
    Defines the differential equations for the coupled spring-mass system.

    Arguments:
        w :  vector of the state variables:
                  w = [x1,y1,x2,y2]
        t :  time
        p :  vector of the parameters:
                  p = [D , gamma,F,Omega, eta, b1, b2]
    """
    x1, y1, x2, y2 = w
    D , gamma,F,Omega, eta, b1, b2 = p

    # Create f = (x1',y1',x2',y2'):
    f = [y1,
         (-b1 * x1 - D * (x1 - x2) -2* gamma * y1-eta*x1**3+F*np.cos(Omega*t)) ,
         y2,
         (-b2 * x2 + D * (x1 - x2)-2* gamma * y2-eta*x2**3) ]
    return f
# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint

# Parameter values

eta = 1.0
F = 1.5
gamma=0.8


Omega = 0.5
D = 1.0

b1 = 0.8
b2 = 0.5

# Initial conditions
# x1 and x2 are the initial displacements; y1 and y2 are the initial velocities
x1 = 0.5
y1 = 0.0
x2 = 2.25
y2 = 0.0

# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 10.0
numpoints = 250

# Create the time samples for the output of the ODE solver.
# I use a large number of points, only because I want to make
# a plot of the solution that looks nice.
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

# Pack up the parameters and initial conditions:
p = [D , gamma,F,Omega, eta, b1, b2]
w0 = [x1, y1, x2, y2]
# Call the ODE solver.
wsol = odeint(vectorfield, w0, t, args=(p,),
              atol=abserr, rtol=relerr)

with open('two_springs.dat', 'w') as f:
    # Print & save the solution.
    for t1, w1 in zip(t, wsol):
        print >> f, t1, w1[0], w1[1], w1[2], w1[3]
# Plot the solution that was generated

from numpy import loadtxt
from pylab import figure, plot, xlabel, grid, hold, legend, title, savefig
from matplotlib.font_manager import FontProperties

t, x1, xy, x2, y2 = loadtxt('two_springs.dat', unpack=True)

figure(1, figsize=(6, 4.5))

xlabel('t')
grid(True)
hold(True)
lw = 1

plot(t, x1, 'b', linewidth=lw)
plot(t, x2, 'g', linewidth=lw)

legend((r'$x_1$', r'$x_2$'), prop=FontProperties(size=16))
title('Mass Displacements for thenCoupled Spring-Mass System')
savefig('two_springs.png', dpi=100)