# Solver error in SciPy/LSODA with a very specific parameter set

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 code I'm using is as follows:

import numpy as np
from pylab import *
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)


This produces an error:

 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.


However, changing "mu" (which is currently 1/70 years, a pretty common value for life expectancy) to either 26550 or 22550 produces no such error, though some values between 25550 and 22550 seem to.

I'm a little puzzled as to why that combination in particular is producing errors. Anyone have some insight?

plot( SIR[:,0] ); show()
I am guessing that the derivative of one of the functions is already machine precision 0 as the error also indicates. Your entire ODE system evolves within tend < 150 and is in steady state during the next following 14850 steps. At some point the derivative between two consecutive time steps is evaluate to 0.