def test_concurrent_ok(self):
        f = lambda t, y: 1.0

        for k in xrange(3):
            for sol in ('vode', 'zvode', 'dopri5', 'dop853'):
                r = ode(f).set_integrator(sol)
                r.set_initial_value(0, 0)

                r2 = ode(f).set_integrator(sol)
                r2.set_initial_value(0, 0)

                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)
                r2.integrate(r2.t + 0.1)

                assert_allclose(r.y, 0.1)
                assert_allclose(r2.y, 0.2)

            for sol in ('dopri5', 'dop853'):
                r = ode(f).set_integrator(sol)
                r.set_initial_value(0, 0)

                r2 = ode(f).set_integrator(sol)
                r2.set_initial_value(0, 0)

                r.integrate(r.t + 0.1)
                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)
                r.integrate(r.t + 0.1)
                r2.integrate(r2.t + 0.1)

                assert_allclose(r.y, 0.3)
                assert_allclose(r2.y, 0.2)

    def __init__(self, lat, lon, h, psi=0, x_earth=0, y_earth=0,
                 integrator='dopri5', use_jac=False, **integrator_params):
        """
        Initialize the equations of the chosen model and selects the
        integrator. Check `scipy.integrate.ode` to see available integrators.

        If use_jac = True the jacobian of the equations is used for the
        integration.
        """
        super().__init__(lat, lon, h, psi, x_earth, y_earth)
        # State vector must be initialized with set_initial_state_vector() method
        self.state_vector = None

        from pyfme.models.euler_flat_earth import lamceq, kaeq, kleq
        self.lamceq = lamceq

        if use_jac:
            from pyfme.models.euler_flat_earth import lamceq_jac, kaeq_jac
            jac_LM_and_AM = lamceq_jac
            jac_att = kaeq_jac
            jac_nav = None  # not implemented
        else:
            jac_LM_and_AM = None
            jac_att = None
            jac_nav = None

        self._ode_lamceq = ode(lamceq, jac=jac_LM_and_AM)
        self._ode_kaqeq = ode(kaeq, jac=jac_att)
        self._ode_kleq = ode(kleq, jac=jac_nav)

        self._ode_lamceq.set_integrator(integrator, **integrator_params)
        self._ode_kaqeq.set_integrator(integrator, **integrator_params)
        self._ode_kleq.set_integrator(integrator, **integrator_params)

def solve(si, y, infile):
    """Conducts the solution step, based on the dopri5 integrator in scipy

    :param si: the simulation info object
    :type si: SimInfo
    :param y: the solution vector
    :type y: np.ndarray
    :param infile: the imported infile module
    :type infile: imported module
    """
    n = ode(f_n).set_integrator('dopri5')
    n.set_initial_value(y0_n(si), si.timer.
                        t0.magnitude)
    n.set_f_params(si)
    th = ode(f_th).set_integrator('dopri5', nsteps=infile.nsteps)
    th.set_initial_value(y0_th(si), si.timer.t0.magnitude)
    th.set_f_params(si)
    while (n.successful()
           and n.t < si.timer.tf.magnitude
           and th.t < si.timer.tf.magnitude):
        si.timer.advance_one_timestep()
        si.db.record_all()
        n.integrate(si.timer.current_time().magnitude)
        update_n(n.t, n.y, si)
        th.integrate(si.timer.current_time().magnitude)
        update_th(th.t, n.y, th.y, si)
    return si.y

    def shooting_function_full(self, r0):
        r"""
        The function used in the shooting algorithm.

        This is the full algorithm from integrating over
        :math:`0 \le \theta \le \pi`. The difference between the
        solution and its derivative at the matching point is the
        error to be minimized.

        Parameters
        ----------

        r0 : list of float
            Initial guess for the horizon radius, as outlined above.

        Returns
        -------

        list of float
            The error at the matching point.
        """
        # First half of the horizon
        H0 = np.array([r0[0], 0.0])
        solver1 = ode(self.expansion)
        solver1.set_integrator("dopri5", atol=1.e-8, rtol=1.e-6)
        solver1.set_initial_value(H0, 0.0)
        solver1.integrate(np.pi / 2.0)
        # Second half of the horizon
        H0 = np.array([r0[1], 0.0])
        solver2 = ode(self.expansion)
        solver2.set_integrator("dopri5", atol=1.e-8, rtol=1.e-6)
        solver2.set_initial_value(H0, np.pi)
        solver2.integrate(np.pi / 2.0)

        return solver1.y - solver2.y

    def __init__(self, integrator="dopri5", model="euler_flat_earth", jac=False, **integrator_params):
        """
        Initialize the equations of the chosen model and selects the
        integrator. Check `scipy.integrate.ode` to see available integrators.

        If jac = True the jacobian of the equations is used for the
        integration.
        """

        if model not in models:
            raise ValueError(
                "The specified model is not available, please \
                             check available models in ..."
            )

        if jac:
            jac_LM_and_AM = euler_flat_earth.jac_linear_and_angular_momentum_eqs
            jac_att = euler_flat_earth.jac_linear_and_angular_momentum_eqs
            jac_nav = None  # not implemented
        else:
            jac_LM_and_AM = None
            jac_att = None
            jac_nav = None

        self._LM_and_AM_eqs = ode(euler_flat_earth.linear_and_angular_momentum_eqs, jac=jac_LM_and_AM)
        self._attitude_eqs = ode(euler_flat_earth.kinematic_angular_eqs, jac=jac_att)
        self._navigation_eqs = ode(euler_flat_earth.navigation_eqs, jac=jac_nav)

        self._LM_and_AM_eqs.set_integrator(integrator, **integrator_params)
        self._attitude_eqs.set_integrator(integrator, **integrator_params)
        self._navigation_eqs.set_integrator(integrator, **integrator_params)

        # State vector must be initialized with set_initial_values() method
        self.state_vector = None

    def s_ode(self,params):
        x0 = np.ones(3)
        s0 = np.zeros(18)
        solver_x = ode(self.dx).set_integrator('dopri5')
        solver_x.set_initial_value(x0,0).set_f_params(params,)

        solver_s = ode(self.ds).set_integrator('dopri5')
        solver_s.set_initial_value(s0,0).set_f_params(params,x0)

        dt = 1
        t1 = 250
        sensitivities = []
        for column in self.cols:
            for param in self.params:
                sensitivities.append('{0},{1}'.format(column,param))
        sol = pandas.DataFrame(np.empty((t1/dt,22)),columns=self.cols+sensitivities)
        i = 0
        #return solver_x,solver_s
        while solver_x.successful() and solver_s.successful() and solver_x.t < t1:
            solver_x.integrate(solver_x.t+dt)
            sol.iloc[i][self.cols] = solver_x.y
            solver_s.set_f_params(params,solver_x.y)
            solver_s.integrate(solver_s.t+dt)
            sol.iloc[i][sensitivities] = solver_s.y
            i += 1
        return sol

def run_trials(z, integrators, times, M):
    # Set up ODE solver
    for integrator in integrators:
        if integrator == 'dop853' or integrator == 'dopri5':
            solver = ode(system)
            solver.set_integrator(integrator, atol=1E-1, rtol=1E-3)
        elif integrator == 'bdf':
            solver = ode(system)
            solver.set_integrator('vode', method=integrator, atol=1E-1, rtol=1E-3, nsteps=2000)

        name = integrator + ' ' + str(times[-1]) + ' ' + str(M)
        try:
            z, c1, c2, t, c1_40km, c2_40km = load_variables(name)
            print "Loaded data for: " + name
        except:
            print "Starting solver: ", integrator, "with times", times
            c1, c2, c1_40km, c2_40km, t = solve(solver, c, times, integrator)
            save_variables(name, z, c1, c2, t, c1_40km, c2_40km)

        # And plot some things
        if times[-1] == 86400.0:
            labels = [str(int(time / 3600.)) + " hours" for time in times[1:]]
            plot_c1(z, c, c1, labels, name)
            plot_c2(z, c, c2, labels, name)
        elif times[-1] == 864000.0:
            plot_40km(t, c1_40km, c2_40km, name)

def _no_collapse_expect_out(num_times,expect_out):
    ##Calculates xpect.values at times tlist if no collapse ops. given
    #  
    #------------------------------------
    opt=odeconfig.options
    if odeconfig.tflag in array([1,10,11]):
        ODE=ode(odeconfig.tdfunc)
        code = compile('ODE.set_f_params('+odeconfig.string+')', '<string>', 'exec')
        exec(code)
    elif odeconfig.tflag==2:
        ODE=ode(_cRHStd)
    elif odeconfig.tflag in array([20,22]):
        ODE=ode(_tdRHStd)
    elif odeconfig.tflag==3:
        ODE=ode(_pyRHSc)
    else:
        ODE = ode(cyq_ode_rhs)
        ODE.set_f_params(odeconfig.h_data, odeconfig.h_ind, odeconfig.h_ptr)
    
    ODE.set_integrator('zvode',method=opt.method,order=opt.order,atol=opt.atol,rtol=opt.rtol,nsteps=opt.nsteps,first_step=opt.first_step,min_step=opt.min_step,max_step=opt.max_step) #initialize ODE solver for RHS
    ODE.set_initial_value(odeconfig.psi0,odeconfig.tlist[0]) #set initial conditions
    for jj in range(odeconfig.e_num):
        expect_out[jj][0]=mc_expect(odeconfig.e_ops_data[jj],odeconfig.e_ops_ind[jj],odeconfig.e_ops_ptr[jj],odeconfig.e_ops_isherm[jj],odeconfig.psi0)
    for k in range(1,num_times):
        ODE.integrate(odeconfig.tlist[k],step=0) #integrate up to tlist[k]
        if ODE.successful():
            state=ODE.y/norm(ODE.y)
            for jj in range(odeconfig.e_num):
                expect_out[jj][k]=mc_expect(odeconfig.e_ops_data[jj],odeconfig.e_ops_ind[jj],odeconfig.e_ops_ptr[jj],odeconfig.e_ops_isherm[jj],state)
        else:
            raise ValueError('Error in ODE solver')
    return expect_out #return times and expectiation values

def _no_collapse_psi_out(num_times,psi_out):
    ##Calculates state vectors at times tlist if no collapse AND no expectation values are given.
    #
    opt=odeconfig.options
    if odeconfig.tflag in array([1,10,11]):
        ODE=ode(odeconfig.tdfunc)
        code = compile('ODE.set_f_params('+odeconfig.string+')', '<string>', 'exec')
        exec(code)
    elif odeconfig.tflag==2:
        ODE=ode(_cRHStd)
    elif odeconfig.tflag in array([20,22]):
        ODE=ode(_tdRHStd)
    elif odeconfig.tflag==3:
        ODE=ode(_pyRHSc)
    else:
        ODE = ode(cyq_ode_rhs)
        ODE.set_f_params(odeconfig.h_data, odeconfig.h_ind, odeconfig.h_ptr)
        
    ODE.set_integrator('zvode',method=opt.method,order=opt.order,atol=opt.atol,rtol=opt.rtol,nsteps=opt.nsteps,first_step=opt.first_step,min_step=opt.min_step,max_step=opt.max_step) #initialize ODE solver for RHS
    ODE.set_initial_value(odeconfig.psi0,odeconfig.tlist[0]) #set initial conditions
    psi_out[0]=Qobj(odeconfig.psi0,odeconfig.psi0_dims,odeconfig.psi0_shape,'ket')
    for k in range(1,num_times):
        ODE.integrate(odeconfig.tlist[k],step=0) #integrate up to tlist[k]
        if ODE.successful():
            psi_out[k]=Qobj(ODE.y/norm(ODE.y,2),odeconfig.psi0_dims,odeconfig.psi0_shape,'ket')
        else:
            raise ValueError('Error in ODE solver')
    return psi_out

def Weightloss_Shooting():
	age, sex =  38. , 'female'
	H, BW    =  1.73, 72.7
	T		 =  12*7.   # Time frame = 3 months    
	
	PALf, EIf = 1.7, 2025
	def EI(t): return EIf
	
	def PAL(t): return PALf
	
	from solution import fat_mass, compute_weight_curve, weight_odes
	F = fat_mass(BW,age,H,sex)
	L = BW-F
	# Fix activity level and determine Energy Intake to achieve Target weight of 145 lbs = 65.90909 kg
	
	init_FL = np.array([F,L])
	EI0, EI1 = 2000, 1950					# Variable Energy Intake
	beta = 65.90909*2.2
	reltol, abstol = 1e-9,1e-8
	dim, iterate = 2, 15
	print '\nj = 1', '\nt = ', EI0
	for j in range(2,iterate):
		print '\nj = ', j, '\nt = ', EI1
		ode_f = lambda t,y: weight_odes(t,y,EI1,PAL(t))
		example1 = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol) 
		example1.set_initial_value(init_FL, 0.) 
		y1 = 2.2*sum(example1.integrate(T) )
		
		if abs(y1-beta)<1e-8: 
			print '\n--Solution computed successfully--'
			break
		# Update
		ode_f = lambda t,y: weight_odes(t,y,EI0,PAL(t))
		example0 = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol) 
		example0.set_initial_value(init_FL, 0.) 
		y0 = 2.2*sum(example0.integrate(T))
		
		EI2 = EI1 - (y1 - beta)*(EI1-EI0)/(y1- y0)
		EI0 = EI1
		EI1 = EI2
		
	# Here we plot the solution 
	ode_f = lambda t,y: weight_odes(t,y,EI1,PAL(t))
	example = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol)
	example.set_initial_value(init_FL, 0.) 
	
	X = np.linspace(0.,T,801)
	Y = np.zeros((len(X),dim))
	Y[0,:] = init_FL
	
	for j in range(1,len(X)): Y[j,:] = example.integrate(X[j]) 
	
	print '\n|y(T) - Target Weight| = ', np.abs(beta-2.2*sum(Y[-1,:]) ),'\n'
	plt.plot(X,2.2*(Y[:,0]+Y[:,1]),'-k',linewidth=2.0)
	plt.axis([-1,T+1,0,170])
	plt.show(); plt.clf()
	return 

 def addMethods(self):
     # override to add _solver function
     ODEsystem.addMethods(self, usePsyco=False)
     if self.haveJacobian():
         self._solver = ode(getattr(self,self.funcspec.spec[1]),
                         getattr(self,self.funcspec.auxfns["Jacobian"][1]))
         self._funcreg['_solver'] = ('self',
              'ode(getattr(self,self.funcspec.spec[1]),' \
                + 'getattr(self,self.funcspec.auxfns["Jacobian"][1]))')
     else:
         self._solver = ode(getattr(self,self.funcspec.spec[1]))
         self._funcreg['_solver'] = ('self', 'ode(getattr(self,' \
                                               + 'self.funcspec.spec[1]))')

def DimensionalForms(BoundConditions, NDParameters, X0=param.X0, d=param.d, rICp=param.rICp, \
                     rhoH=param.rhoH, alpha=param.alpha, rhoD=param.rhoD, g=param.g, hminus=param.hmin, \
                     hmaxi=param.hmax, dh=param.dt, geom="cart", AbsTole=param.AbsTol, RelTole=param.AbsTol):
    """ Integration RK (de type Dormand Prince 4-5) of the equations system for an horizontal F-layer.
    resolutionSystem(y0, t0, tmax, dt, *arg)
    Default parameters for the others arguments are :
    PeT=1.e6,
    MX=0.17,
    MP=0.016,
    Veff=1.e7,
    gamma=7.2,
    Atol=1.e-10,
    Rtol=1.e-12
    The result is given in fully dimensional form. 
    """

    if geom == "cart":
        sol = ode(systemdiff.cart).set_integrator(
            'dopri5', atol=param.AbsTol, rtol=param.AbsTol, nsteps=1e7)
    elif geom == "spher":
        sol = ode(systemdiff.spher).set_integrator(
            'dopri5', atol=param.AbsTol, rtol=param.AbsTol, nsteps=1e7)

    sol.set_initial_value(BoundConditions, param.hmin)
    sol.set_f_params(NDParameters)  # (PeT, Mx, Mp, Veff, gamma)
    X, dXdh, phiVs,  h = BoundConditions[0] * np.ones(1), BoundConditions[1] * np.ones(1), \
      BoundConditions[2] * np.ones(1), param.hmin * np.ones(1)

    if param.dt > 0:
        while sol.successful() and sol.t < param.hmax:
            sol.integrate(sol.t + param.dt)
        # print("%g %g" % (r.t, r.y))
            X = np.append(X, sol.y[0])
            dXdh = np.append(dXdh, sol.y[1])
            phiVs = np.append(phiVs, sol.y[2])
            h = np.append(h, sol.t)
    else:
        while sol.successful() and sol.t > param.hmax:
            sol.integrate(sol.t + param.dt)
        # print("%g %g" % (r.t, r.y))
            X = np.append(X, sol.y[0])
            dXdh = np.append(dXdh, sol.y[1])
            phiVs = np.append(phiVs, sol.y[2])
            h = np.append(h, sol.t)

    # back to fully dimensional:
    x = X0 * X
    z = h * d
    T = param.Tliq0 * (1 - NDParameters['Mp'] * h - X * NDParameters['Mx'])

    return z, x, phiVs, T

def run():
    data_x = []
    data_y = []

    dt = 0.1
    tmax = 10
    state = None

    m = model.model_push()
    ir = sciint.ode(m.step).set_integrator('vode', atol=0.01)
    ir.set_initial_value(m.init, 0)

    while state is None or state[2] < m.x_h and ir.t < tmax:
        state = ir.integrate(ir.t + dt/1000)
        data_x.append(ir.t)
        data_y.append(m.M)

    if ir.t >= tmax:
        print('front failed')
        plt.ylim(0, m.Mmax + 10)
        plt.plot(data_x, data_y)
        plt.xlabel('Время t [с]')
        plt.ylabel('Момент М(t) [Н м]')
        plt.show()
        return None

    t1 = ir.t
    m = model.model_torque()
    ir = sciint.ode(m.step).set_integrator('vode', atol=0.01)
    m.init = [
        state[0], 0,
        state[2] - m.L * cos(asin(((state[3] - m.y_f)/m.L))), m.r_k + m.y_f
    ]
    state2 = state[:]
    ir.set_initial_value(m.init, 0)
    state = None
    while state is None or state[2] < m.x_h and ir.t + t1 < tmax:
        state = ir.integrate(ir.t + dt/1000)
        data_x.append(ir.t + t1)
        data_y.append(m.M)

    if ir.t + t1 >= tmax:
        print('rear failed')
        return None
    print('all pass')
    plt.ylim(0, m.Mmax + 10)
    plt.plot(data_x, data_y)
    plt.xlabel('Время t [с]')
    plt.ylabel('Момент М(t) [Н м]')
    plt.show()
    return (m.init, state2)

    def test_concurrent_fail(self):
        for sol in ('vode', 'zvode'):
            f = lambda t, y: 1.0

            r = ode(f).set_integrator(sol)
            r.set_initial_value(0, 0)

            r2 = ode(f).set_integrator(sol)
            r2.set_initial_value(0, 0)

            r.integrate(r.t + 0.1)
            r2.integrate(r2.t + 0.1)

            assert_raises(RuntimeError, r.integrate, r.t + 0.1)

def Cannon_Shooting(t0,t1,za=0.,va=45,nu=.0003): 
	a, b = 0., 195
	beta =  0.
	
	# t0,t1 = np.pi/6., np.pi/7
	dim, iterate = 3,40
	reltol, abstol = 1e-9,1e-8
	
	# Initial_Conditions = np.array([z=za,v=va,phi=t])
	g = 9.8067
	def ode_f(x,y): 
		# y = [z,v,phi]
		return np.array([np.tan(y[2]), -(g*np.sin(y[2]) + nu*y[1]**2.)/(y[1]*np.cos(y[2])), 
		-g/y[1]**2.])
	
	
	print '\nj = 1'
	print 't = ', t0
	for j in range(2,iterate):
		print '\nj = ', j
		print 't = ', t1
		example1 = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol) 
		example1.set_initial_value(np.array([za,va,t1]),a) 
		y1 = example1.integrate(b)[0]
		
		if abs(y1-beta)<1e-8: 
			print '\n--Solution y computed successfully--'
			break
		# Update
		example0 = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol) 
		example0.set_initial_value(np.array([za,va,t0]),a) 
		y0 = example0.integrate(b)[0]
		
		t2 = t1 - (y1 - beta)*(t1-t0)/(y1- y0)
		t0 = t1
		t1 = t2
		
	# Here we plot the solution 
	example = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol)
	example.set_initial_value(np.array([za,va,t1]),a) 
	X = np.linspace(a,b,801)
	Y = np.zeros((len(X),dim))
	Y[0,:] = np.array([0.,0.5,t1])
	
	for j in range(1,len(X)): Y[j,:] = example.integrate(X[j]) 
	
	print '\n|y(b) - beta| = ', np.abs(beta-Y[-1,0]),'\n'
	# plt.plot(X,Y[:,0],'-k')
	return X,Y

def Exercise1(): 
# y'' +4y = -9sin(x), y(0) = 1., y(3*pi/4.) = -(1.+3*sqrt(2))/2., y'(0) = -2
# Exact Solution: y(x) = cos(2x) + (1/2)sin(2x) - 3sin(x)
	
	a, b = 0., 3*np.pi/4.
	alpha, beta =  1., -(1.+3*np.sqrt(2))/2.
	t0,t1 = 3., 1.
	dim, iterate = 2,10
	reltol, abstol = 1e-9,1e-8
	def ode_f(x,y): return np.array([y[1] , -4.*y[0] - 9.*np.sin(x)])
	
	
	print '\nj = 1','\nt = ', t0
	
	for j in range(2,iterate):
		print '\nj = ', j,'\nt = ', t1
		
		if j >2: 
			y0 = y1			# Update
		else:
			example = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol) 
			example.set_initial_value(np.array([alpha,t0]),a) 
			y0 = example.integrate(b)[0]
			
		example = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol) 
		example.set_initial_value(np.array([alpha,t1]),a) 
		y1 = example.integrate(b)[0]
		if abs(y1-beta)<1e-8: 
			print '\n--Solution y computed successfully--',\
			'\n|y(b) - beta| = ',np.abs(beta-y1),'\n'
			break
		
		t2 = t1 - (y1 - beta)*(t1-t0)/(y1- y0)
		t0 = t1
		t1 = t2
	
	# Plots the solution
	example = ode(ode_f).set_integrator('dopri5',atol=abstol,rtol=reltol)
	example.set_initial_value(np.array([alpha,t1]),a) 
	X = np.linspace(a,b,401)
	Y = np.zeros((len(X),dim))
	Y[0,:] = np.array([alpha, t1])
	
	for j in range(1,len(X)): 
		Y[j,:] = example.integrate(X[j]) 
	plt.plot(X,Y[:,0],'-k')
	plt.show()
	return 

def implicit_black_box(propensities, V, X, w, h, deter_vector, stoc_positions, positions, valid, deriv):

    # Adjustment for systems reaching steady state

    temp = derivative_G(propensities, V, X, w, deter_vector, stoc_positions, positions, valid)
    # pdb.set_trace()
    valid_adjust_pos = np.where(np.sum(np.abs(temp), axis=0) < 1e-10, True, False)

    valid_adjust = valid[:, :]
    valid_adjust[valid_adjust_pos, :] = False

    # print(" Reached Steady State %d"%(np.sum(valid_adjust_pos)))

    from scipy.integrate import ode

    # pdb.set_trace()
    deter_ode = ode(f).set_integrator("lsoda", method="adams", with_jacobian=False)
    deter_ode.set_initial_value(X[deter_vector, :].flatten(), 0).set_f_params(
        [propensities, V, X, deter_vector, stoc_positions, positions, valid_adjust, w]
    )

    # pdb.set_trace()
    while deter_ode.successful() and deter_ode.t < h:
        deter_ode.integrate(h)

        # print("Black Box: \n"+ str(deter_ode.y))

        # print("iterator : \n:"+str(next_X[deter_vector,:]))
    X[deter_vector, :] = deter_ode.y.reshape((np.sum(deter_vector), X.shape[1]))

    # Another adjust to compensate for non negative
    X = np.where(X < 0.0, 0.0, X)

    return X

    def numerical_integrator(self, t0, tf, dt, initial_conditions, solver, verbose):
        """
        the numerical_integrator() method integrates the ODE of the dynamic system using a numerical solver
        """
        f = self._ode_RHS
        if initial_conditions:
            y0 = initial_conditions
        else:
            y0 = self.initial_conditions

        MdFBA_ode = ode(f).set_integrator(solver)
        MdFBA_ode.set_initial_value(y0, t0)

        t = [t0]
        y = [y0]

        while MdFBA_ode.successful() and MdFBA_ode.t < tf:
            MdFBA_ode.integrate(MdFBA_ode.t + dt)

            t.append(MdFBA_ode.t)
            y.append(MdFBA_ode.y)

            if verbose:
                print MdFBA_ode.t

        t = numpy.array(t)
        y = numpy.array(y)
        return t, y

def run_ODE(f, a_in, C, D, num_of_variables, T = 10, dt = 0.01, y0 = None):
    '''Run the ODE for the given set of equations and record the outputs.

    Args:
        f (function): Evolution of the system.

        a_in (function): inputs as a function of time.

        C,D (matrices): matrices to use to obtain output from system state and
            input.

        num_of_variables (int): number of variables the system has.

        T (optional[positive float]): length of simulation.

        dt (optional[float]): time step used by the simulation.

    Returns:
        (array): 
        An array Y of outputs.

    '''
    if y0 is None:
        y0 = np.asmatrix([0.]*num_of_variables).T

    r = ode(f).set_integrator('zvode', method='bdf')
    r.set_initial_value(y0, 0.)
    Y = []
    while r.successful() and r.t < T:
        Y.append(C*r.y+D*a_in(r.t))
        r.integrate(r.t+dt)
    return Y

    def run(self, start, velocity):
        # state data structures
        self.t=zeros(self.num_steps)
        self.Y=zeros([self.num_steps, 4])
        self.Y[0] = array([start[0],start[1],velocity[0],velocity[1]])

        r = ode(self.deriv)
        r.set_integrator('vode', rtol=1e-9, atol=1e-6)
        r.set_initial_value(self.Y[0],t=0.0)

        time_idx=1
        while r.successful() and r.t < self.max_time:
            r.integrate(time_idx*self.time_step)
            self.t[time_idx] = r.t
            self.Y[time_idx] = r.y
            time_idx = time_idx + 1

            if self.debug:
                print ("time=%10.3f: x=%g y=%g vx=%g vy=%g" % \
                      (r.t, r.y[0], r.y[1], r.y[2], r.y[3]))

            if r.y[1]<0:
                break
            
        # simulation done: truncate arrays 
        self.Y.resize([time_idx, self.Y.shape[1]])
        self.t.resize([time_idx])