[SciPy-User] Odd scipy.integrate error at a specific parameter value

[Lofgren, Eric]

Cross posted to http://scicomp.stackexchange.com

I'm implementing a very simple Susceptible-Infected-Recovered model with a steady population for an idle side project - normally a pretty trivial task. But I'm running into solver errors using either PysCeS or SciPy, both of which use lsoda as their underlying solver. This only happens for particular values of a parameter, and I'm stumped as to why. The system of differential equations is as follows:

dS/dt = -beta*S*I/N + mu*N - mu*S
dI/dt = beta*S*I/N - gamma*I - mu*I
dR/dt = gamma*I - mu*R

Where N = S + I + R.

The code I'm using for this is this:

####
import numpy as np
import scipy.integrate as spi

#Parameter Values
S0 = 99.
I0 = 1.
R0 = 0.
N0 = S0 + I0 + R0
PopIn= (S0, I0, R0, N0)
beta= 0.50     
gamma=1/10.  
mu = 1/25550.
t_end = 15000.
t_start = 1.
t_step = 1.
t_interval = np.arange(t_start, t_end, t_step)

#Solving the differential equation. Solves over t for initial conditions PopIn
def eq_system(PopIn,t):
    '''Defining SIR System of Equations'''
    #Creating an array of equations
    Eqs= np.zeros((3))
    Eqs[0]= -beta * (PopIn[0]*PopIn[1]/PopIn[3]) - mu*PopIn[0] + mu*PopIn[3]
    Eqs[1]= (beta * (PopIn[0]*PopIn[1]/PopIn[3]) - gamma*PopIn[1] - mu*PopIn[1])
    Eqs[2]= gamma*PopIn[1] - mu*PopIn[2]
    return Eqs

SIR = spi.odeint(eq_system, PopIn, t_interval)
####

The problem is that using 1/25550, which is ~ 1/70 years in days, a pretty common value for a birth/death rate, the solver throws an error at around t = 7732, producing the following:

 lsoda--  at current t (=r1), mxstep (=i1) steps   
       taken on this call before reaching tout     
      In above message,  I1 =       500
      In above message,  R1 =  0.7732042715460E+04
Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.

This is true for values between about 22550 and 25550, but setting mu = 1/26550 or 1/22550 ends up producing no such error, and indeed showing exactly the kind of behavior one would expect.

Does anyone have an idea of what's causing scipy.integrate/lsoda to grind to a halt with that parameter value?

Thanks,

Eric