<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 SIR Model
The example describes how, given a contagious illness, the population varies among the three groups of the indeividuals 
that are susceptible to the illness ($S$), infected by the illness ($I$) and recovered from the illness ($R$).
Nonlinear Ordinary Differential Equations
for 
$$ S\overset{\beta I}{\longrightarrow}I\overset{\gamma}{\longrightarrow}R$$
\begin{eqnarray}
  \frac {d S}{d t} &=& -\beta S I\
   \frac {d I}{d t} &=& \beta S I -\gamma I\
   \frac {d R}{d t} &=& \gamma I
\end{eqnarray}
Step1:  Solve system of ODE's by Scipy
Procedures
define time partition, t_range, initials, INPUT
Step2:  Reproduction Ratio
<ol>
<li> $\frac{d I}{d t}$ 
Step3:  Note
Facts
Step4:  Conclusion
After $3.0758505917569066$ units of time, it  would be under control.
Computer Practice
Make an interactive demo with input, $\beta,\gamma,R_0,I_0$, and $T$ output.
Zombie Apocalypse
Philip Munz, Ioan Hudea, Joe Imad and Robert J. Smith 2009
\begin{eqnarray}
  \frac {d S}{d t} &=& -\beta S Z\
   \frac {d Z}{d t} &=& (\beta-\alpha) S Z +\zeta R\
   \frac {d R}{d t} &=& \alpha S Z -\zeta R\
\end{eqnarray}
Encounters between zombies and live humans can result in the human becoming a zombie, $\beta S Z$ ;
the zombie becoming dead, $\alpha S Z$;
The dead sometimes spontaneously turn into zombies, $\zeta R$.
Step5:  Ebola
Step6:  Governed Differential Equations
\begin{eqnarray}
  \frac {d S}{d t} &=& -\frac{\beta S I}{N}\
  \frac {d E}{d t} &=&\frac{\beta S I}{N}- k E\
   \frac {d I}{d t} &=& k E -\gamma I\
   \frac {d R}{d t} &=& \gamma I \
   \frac {d C}{d t} &=& k E\
\end{eqnarray}
where $C(t)$ is not an epidemiological state but serves to keep track of the cumulative number of Ebola cases from the time of onset of symptoms.
$1/k$ Mean incubation period (days)
$1/\gamma$ Mean infectious period (days)
$\beta=\frac{\beta_0+\beta_1}{2}$
Case for Congo 1995
N=564500
$\beta_0=0.33,\beta_1=0.09$
$1/k=5.3,1/\gamma=5.61$
Step7:  Intervention Case
After diease breakout, there ahould be some intervention strategies made to control the spread of Ebola.
The transmission rate, $\beta$, could be modeled reasonably as follows
Step8:  Symboilc Solution dor ODE's 
Up to now, we only solve ODE's by numeric schemes (approximately); can we solve them theoretically? Just as well-know fact
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numpy import log
from scipy import integrate
from scipy.optimize import fsolve
import scipy.linalg as la
import scipy.integrate as spi
from IPython.html.widgets import interact, interactive, fixed
from IPython.html import widgets
#from matplotlib import rc
#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
## for Palatino and other serif fonts use:
#rc('font',**{'family':'serif','serif':['Palatino']})
#rc('text', usetex=True)
S0=0.995;I0=.005;R0=0.
beta=0.;gamma=0.35
n=1000
T=20
dt=T/float(n)
t=np.linspace(0,T,n+1)
plt.xlim(0,20)
plt.ylim(0,1.1)
beta=0.75
S,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
S[0]=S0;I[0]=I0;R[0]=R0
for i in range(len(t)-1):
    S[i+1] = S[i]-dt*beta*S[i]*I[i]
    I[i+1] = I[i]+dt*(beta*S[i]*I[i]-gamma*I[i]) 
    #R[i+1] = 1.-S[i+1]-I[i+1]
R=1.-S-I
plt.plot(t,S,t,I,t,R) 
plt.legend(['S','I','R'])
plt.grid(True)
plt.title(r'$\beta=0.75,\gamma=0.35$',size=20)
    
Explanation: SIR Model
The example describes how, given a contagious illness, the population varies among the three groups of the indeividuals 
that are susceptible to the illness ($S$), infected by the illness ($I$) and recovered from the illness ($R$).
Nonlinear Ordinary Differential Equations
for 
$$ S\overset{\beta I}{\longrightarrow}I\overset{\gamma}{\longrightarrow}R$$
\begin{eqnarray}
  \frac {d S}{d t} &=& -\beta S I\
   \frac {d I}{d t} &=& \beta S I -\gamma I\
   \frac {d R}{d t} &=& \gamma I
\end{eqnarray}
End of explanation
plt.xlim(0,20)
plt.ylim(0,1.1)
beta=0.75
S,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
S[0]=S0;I[0]=I0;R[0]=R0
for i in range(len(t)-1):
    S[i+1] = S[i]-dt*beta*S[i]*I[i]
    I[i+1] = I[i]+dt*(beta*S[i]*I[i]-gamma*I[i]) 
    #R[i+1] = 1.-S[i+1]-I[i+1]
R=1.-S-I
def diff_eqs(INP,t):  
    '''The main set of equations'''
    Y=np.zeros((2))
    V = INP    
    Y[0] = - beta * V[0] * V[1]
    Y[1] = beta * V[0] * V[1] - gamma * V[1] 
    return Y   # For odeint
t_start = 0.0; t_end = T; t_inc = dt
t_range = np.arange(t_start, t_end+t_inc, t_inc)
INPUT=(S0,I0)
RES = spi.odeint(diff_eqs,INPUT,t_range)
Rec=1. - (RES[:,0]+RES[:,1])
plt.plot(t,RES[:,0],t,RES[:,1],t,Rec) 
plt.legend(['S','I','R'])
plt.grid(True)
plt.title(r'$\beta=0.75,\gamma=0.35$',size=20)
    
beta=2.
def vis2(beta,gamma):
    fig, ax = plt.subplots(figsize=(8,4),
                           subplot_kw={'axisbg':'#EEEEEE',
                                       'axisbelow':True})
    S0=0.995;I0=.005;R0=0.
    n=1000
    T=20
    dt=T/float(n)
    t=np.linspace(0,T,n+1)
    
    ax.set_xlim(0,T)
    ax.set_ylim(0,1.1)
    S,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
    S[0]=S0;I[0]=I0;R[0]=R0
    for i in range(len(t)-1):
        S[i+1] = S[i]-dt*beta*S[i]*I[i]
        I[i+1] = I[i]+dt*(beta*S[i]*I[i]-gamma*I[i]) 
        R[i+1] = 1.-S[i+1]-I[i+1]
    ax.plot(t,S)
    ax.plot(t,I)
    ax.plot(t,R)
    ax.set_title("SIR model")
    ax.set_xlabel("Time (days)")
    ax.set_ylabel("Percent of Population")
    ax.legend(('S','I','R'))
    ax.grid(True)
    return fig
i = interact(vis2,
         beta=widgets.FloatSliderWidget(min=0.0, max=5.0, step=0.1, value=2.0, description="beta"),
         gamma=widgets.FloatSliderWidget(min=1.0, max=5, step=0.1, value=1.5, description="gamma"),
         )
beta=2.
gamma=1.5
S0=0.95
I0=0.05
R0=0.
def vis2(beta,gamma,S0,I0):
    fig, ax = plt.subplots(figsize=(8,4),
                           subplot_kw={'axisbg':'#EEEEEE',
                                       'axisbelow':True})
    n=1000
    T=40
    R0=1-S0-I0
    dt=T/float(n)
    
    t=np.linspace(0,T,n+1)
    R0=gamma/beta
    ax.set_xlim(0,T)
    ax.set_ylim(0,1.1)
    S,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
    S[0]=S0;I[0]=I0;R[0]=R0
    for i in range(len(t)-1):
        S[i+1] = S[i]-dt*beta*S[i]*I[i]
        I[i+1] = I[i]+dt*(beta*S[i]*I[i]-gamma*I[i]) 
        R[i+1] = 1.-S[i+1]-I[i+1]
    ax.plot(t,S)
    ax.plot(t,I)
    ax.plot(t,R)
    ax.set_title('SIR model $R_0=$%f' %R0)
    ax.set_xlabel("Time (days)")
    ax.set_ylabel("Percent of Population")
    ax.legend(('S','I','R'))
    ax.grid(True)
    return fig
i = interact(vis2,
         beta=widgets.FloatSliderWidget(min=0.0, max=5.0, step=0.1, value=2.0, description="beta"),
         gamma=widgets.FloatSliderWidget(min=1.0, max=2.0, step=0.1, value=1.5, description="gamma"),
         S0=widgets.FloatSliderWidget(min=0.9, max=1.0, step=0.01, value=0.95, description="S0"),
         I0=widgets.FloatSliderWidget(min=0., max=0.1, step=0.01, value=0.01, description="I0"),
         )
Explanation: Solve system of ODE's by Scipy
Procedures
define time partition, t_range, initials, INPUT:(S0,I0);
define Governed equation,diff_eqs;
solved by scipy.integrate.odeint(diff_eqs, INPUT, t_range)
End of explanation
beta=2.
gamma=1.5
S0=0.3
I0=0.05
R0=0.
def SI3(beta,gamma):
    fig, ax = plt.subplots(figsize=(8,4),
                           subplot_kw={'axisbg':'#EEEEEE',
                                       'axisbelow':True})
    n=1000
    St=np.linspace(S0,1,n+1)
    
    ax.set_xlim(0,1)
    ax.set_ylim(0,1.1)
    I=gamma*log(St)/beta-St+I0+S0-gamma*log(S0)/beta
    ax.plot(St,I)
    ax.set_title("S-I Relation")
    ax.set_xlabel("S")
    ax.set_ylabel("I")
    ax.grid(True)
    return fig
i = interact(SI3,
         beta=widgets.FloatSliderWidget(min=0.0, max=10.0, step=0.1, value=2.0, description="beta"),
         gamma=widgets.FloatSliderWidget(min=1.0, max=10., step=0.1, value=1.5, description="gamma"),
         )
Explanation: Reproduction Ratio
<ol>
<li> $\frac{d I}{d t}$ : rate of changes (growth velocity) of infectives. 
<li> The number of infectives grows as long as the proportion of susceptibles, $S(t)$, 
      in the population exceeds $\frac{\gamma}{\beta}$, i.e. 
$$ \frac{d I}{d t}=\beta S I -\gamma I\ge(\beta\frac{\gamma}{\beta}-\gamma)\cdot I\ge0$$
</ol>
Defintion
$R_0=\frac{\gamma}{\beta}$ is called repoduction ratio, and
$R_0>1$: epidemic is out of control.
$R_0<1$: epidemic is under control.
What relation between $S$ and $I$
Theoretically, it's hard to solve from $S-I-R$ ODE's.
Fortunately, we can reduce the $S-I$ DE's into:
$$ \frac { d I}{d S}=\frac{\frac{d I}{d S}}{\frac{d S}{d I}}=\frac{\beta I S-\gamma I}{-\beta I S}
    =-1+\frac{\gamma}{\beta S}$$
And 
$$ I=\int\left(\frac{\gamma}{\beta S} -1\right) d S=\frac{\gamma}{\beta} \ln S -S+C$$
Replcing the initial values, $S(0),I(0)$, gets the constant $C$:
$$ C=I(0)-\frac{\gamma}{\beta}\ln(S(0))+S(0)$$
End of explanation
b=8.
g=2.3
S0=0.99
I0=0.01
R0=0.
C0= I0-g/b*log(S0)+S0
print C0
def func(x):
        return -g/b*np.log(x)+x+I0-C0
SS=fsolve(func,0.01);print(SS)
def xv0(x):
    return 0*x
def intf(x):
    return 1/(b*x*x-b*C0*x-g*x*log(x))
xv=np.linspace(SS,S0,100)
xv2=np.linspace(0,1,101)
plt.ylim(-50,5)
plt.text(0.2,-20,r'$\int^{S_1}_{S_0}\frac{ d S }{-\gamma S \ln S +\beta S^2-\beta C S}$',size=20)
plt.text(0.2,-30,r'$S_0=0.99,S_1=0.0358$',size=12)
plt.plot(xv2,xv0(xv2))
plt.plot(xv,intf(xv),'b--')
plt.fill_between(xv, intf(xv),0, facecolor='red', alpha=0.1)
x2 = lambda v: 1/(b*v*v-b*C0*v-g*v*log(v))
integrate.quad(x2,S0,SS)
Explanation: Note
Facts:
$I=0\to S,I,R$ are contstants: $I$ is an equivalent state, i.e. nothing will happen.
If case starts out at $I(0)>0$, $\frac{ d S}{d t}<0$ implies $S$ being decreasing. 
$ I=\frac{\gamma}{\beta}\ln S -S +C \to -\infty$ as $S\to0$.
Conclusion
Epidemic will end as $I\to 0$ while $S$ approaches some positive term.
While $I$ is at equilibium, $S(t)=S_T$ for some $t\ge t_S$
Questions
Observe the relation between $S$ and $I$ by making plots of $S-I$
Problem: How long is $t_S$?
To get more precisely, i.e. when would it be in control, duration time for approaching $S(t)$ has to to be estimested. 
Replacing $I=I(S)$ above gets:
\begin{eqnarray}
 \frac{d S}{d t} &= &-\beta I S \cr
  &= &-\gamma S \ln S +\beta S^2-\beta C S 
\end{eqnarray} 
This implies
$$\int^{S_T}_{S_0}\frac{ d S }{-\gamma S \ln S +\beta S^2-\beta C S} =\int^T_0 d t =T$$
Solution Processes
It's true: epidemic will explose and then enventually be under control.
It is more interesting to know when  $I$ come back to $I_0$ again than expect to get $I=0$.
As we note above, when to get $S=S_1$ while $I=I_0$ come back again?
Answer
This is the key point:
\begin{eqnarray}
C&=&I(t)-\frac{\gamma}{\beta}\ln S(t)+S(t) \
  &=& I(0)-\frac{\gamma}{\beta}\ln(S(0))+S(0) \
  &=& I(0)-\frac{\gamma}{\beta}\ln S_1 +S_1
\end{eqnarray}
Replacing the upper limit in the integral and get the time, $T$, when the epidemic will be under control:
\begin{eqnarray}
 T&=&\int^T_0 d t \
  &=&\int^{S_1}_{S_0}\frac{ d S }{-\gamma S \ln S +\beta S^2-\beta C S}\qquad(\text{ where } S_1<S_0)
\end{eqnarray}
End of explanation
beta=1.
alpha=4.
zeta=0.2
S0=.99
Z0=0.01
R0=0.0
T=100.
n=10000
t=np.linspace(0,T,n+1)
dt=t[1]-t[0]
S,Z,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
S[0]=S0;Z[0]=Z0;R[0]=R0
for i in range(len(t)-1):
    S[i+1] = S[i]-dt*beta*S[i]*Z[i]
    Z[i+1] = Z[i]+dt*((beta-alpha)*S[i]*Z[i]+zeta*R[i]) 
    R[i+1] = 1.-S[i+1]-Z[i+1]
plt.plot(t,S,t,Z,t,R)
plt.legend(('S','Z','R'))
print('S[100]','Z[100]','R[100]:',S[-1],Z[-1],R[-1])
t[:5]
Explanation: Conclusion
After $3.0758505917569066$ units of time, it  would be under control.
Computer Practice
Make an interactive demo with input, $\beta,\gamma,R_0,I_0$, and $T$ output.
Zombie Apocalypse
Philip Munz, Ioan Hudea, Joe Imad and Robert J. Smith 2009
\begin{eqnarray}
  \frac {d S}{d t} &=& -\beta S Z\
   \frac {d Z}{d t} &=& (\beta-\alpha) S Z +\zeta R\
   \frac {d R}{d t} &=& \alpha S Z -\zeta R\
\end{eqnarray}
Encounters between zombies and live humans can result in the human becoming a zombie, $\beta S Z$ ;
the zombie becoming dead, $\alpha S Z$;
The dead sometimes spontaneously turn into zombies, $\zeta R$.
End of explanation
mu=1/(70*365.0)
beta=520/365.0
k=1/14.0
gamma=1/7.0
ND=10*365.0
TS=1.0
S0=0.1
E0=1e-4
I0=1e-4
INPUT = (S0, E0, I0)
def diff_eqs(INP,t):  
	'''The main set of equations'''
	Y=np.zeros((3))
	V = INP    
	Y[0] = mu - beta * V[0] * V[2] - mu * V[0]
	Y[1] = beta * V[0] * V[2] - k * V[1] - mu * V[1]
	Y[2] = k * V[1] - gamma * V[2] - mu * V[2]
	return Y   # For odeint
t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range)
Rec=1. - (RES[:,0]+RES[:,1]+RES[:,2])
#print RES
#Ploting
plt.subplot(311)
plt.plot(RES[:,0], '-g', label='Susceptibles')
plt.title('SEIR with Births and Deaths')
plt.xlabel('Time')
plt.ylabel('Susceptibles')
plt.subplot(312)
plt.plot(RES[:,1], '-m', label='Exposed')
plt.plot(RES[:,2], '-r', label='Infectious')
plt.legend(loc=0)
plt.xlabel('Time')
plt.ylabel('Infected')
plt.subplot(313)
plt.plot(Rec, '-k', label='Recovereds')
plt.xlabel('Time')
plt.ylabel('Recovereds')
beta=520/365.
gamma=1/7.0
#mu=1/(70*365.0)
mudays=70
def vis3(mudays,beta,gamma):
    #mu=1/(70*365.0)
    #beta=520/365.0
    k=1/14.0
    #gamma=1/7.0
    mu=1/float(mudays*365)
    ND=10*365.0
    TS=1.0
    S0=0.1
    E0=1e-4
    I0=1e-4
    INPUT = (S0, E0, I0)
    def diff_eqs(INP,t):  
        '''The main set of equations'''
        Y=np.zeros((3))
        V = INP    
        Y[0] = mu - beta * V[0] * V[2] - mu * V[0]
        Y[1] = beta * V[0] * V[2] - k * V[1] - mu * V[1]
        Y[2] = k * V[1] - gamma * V[2] - mu * V[2]
        return Y   # For odeint
    
    t_start = 0.0; t_end = ND; t_inc = TS
    t_range = np.arange(t_start, t_end+t_inc, t_inc)
    RES = spi.odeint(diff_eqs,INPUT,t_range)
    Rec=1. - (RES[:,0]+RES[:,1]+RES[:,2])
    
    fig, ax = plt.subplots()
    plt.subplot(311)
    plt.plot(RES[:,0], '-g', label='Susceptibles')
    plt.title('SEIR with Births and Deaths,mu=1/mudys/365')
    plt.xlabel('Time')
    plt.ylabel('Susceptibles')
    plt.subplot(312)
    plt.plot(RES[:,1], '-m', label='Exposed')
    plt.plot(RES[:,2], '-r', label='Infectious')
    plt.legend(loc=0)
    plt.xlabel('Time')
    plt.ylabel('Infected')
    plt.subplot(313)
    plt.plot(Rec, '-k', label='Recovereds')
    plt.xlabel('Time')
    plt.ylabel('Recovereds')
    
    return ax
i = interact(vis3,
         mudays  =widgets.FloatSliderWidget(min=1, max=365, step=1, value=70, description="mudays"),    
         beta=widgets.FloatSliderWidget(min=0.0, max=5.0, step=0.1, value=520/365., description="beta"),
         gamma=widgets.FloatSliderWidget(min=0.0, max=5, step=0.1, value=1/7., description="gamma"),
         )
Explanation: Ebola: SEIR Model
Ref:
G.Chowell, N.W.Hengartner,C.Castillo-Chavez,P.W.Fenimore,J.M.Hyman, The basic reproductive number of Ebola and the effects of public health measures: the cases of Congo and Uganda,  Journal of Theoretical Biology, 2004
Anderson, R.M., May, R.M., 1991. Infectious Diseases of Humans. Oxford University Press, Oxford.
Brauer, F., Castillo-Chavez, C., 2000. Mathematical Models in Sartwell, P.E., 1966. The distribution of incubation periods and the PopulationBiologyandEpidemiology.Springer,NewYork.
Matt J. Keeling & Pejman Rohani,
  <a href="http://wiki.deductivethinking.com">Modeling Infectious Diseases in Humans and Animals</a>
Susceptible individuals in class $S$ in contact with the
virus enter the exposed class $E$ at the per-capita rate $\beta I/N$; 
where $\beta$ is transmission rate per person per day, 
$N$is the total effective population size,
and $I/N$ is the probability that a contact is made with a infectious individual (i.e. uniform mixing is assumed). 
Exposed individuals undergo an average incubation period (assumed asymptomatic and uninfectious) 
of $1/k$ days before progressing to the infectious class $I$: Infectious individuals move to the $R$-class
(death or recovered) at the per-capita rate $\gamma$, as figure:
<img src="images/SEIR2.png">
<img src="images/SEIR1.png">
Simple SEIR Simulation ( with births and deaths)
\begin{eqnarray}
\frac{dS}{dt} &= \mu-(\beta I+\mu) S \
\frac{dE}{dt} &= \beta S I-(\mu+k) E \
\frac{dI}{dt} &= k E-(\mu+\gamma) I \
\frac{dR}{dt} &= \gamma I-\mu E 
\end{eqnarray}
End of explanation
N=64500.
S0=0.9*N;E0=0.0;I0=N-S0-E0;R0=0.;C0=0.;
beta0=0.33;beta1=0.09;
beta=beta0;gamma=1/5.61
k=1/5.3
n=10000
T=100
dt=T/float(n)
t=np.linspace(0,T,n+1)
S,E,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
S[0]=S0;E[0]=E0;I[0]=I0;R[0]=R0
for i in range(len(t)-1):
    S[i+1] = S[i]-dt*beta*S[i]*I[i]/N
    E[i+1] = E[i]+dt*(beta*S[i]*I[i]/N-k*E[i])
    I[i+1] = I[i]+dt*(k*E[i]-gamma*I[i]) 
    R[i+1] = N-S[i+1]-I[i+1]-E[i+1]
plt.plot(t,S/N,t,I/N,'r--',t,E/N,'y',t,R/N,'k-') 
#plt.ylim(0,1)
plt.legend(['Suspective','Infected','Exposed','Recoverd'])
beta/gamma
J=np.matrix([[-k, beta],[k, -gamma]]);J
import scipy.linalg as la
r=np.max(np.real(la.eigvals(J)));r
#r=0.07
1+(r*r+(k+gamma)*r)/k/gamma
Explanation: Governed Differential Equations
\begin{eqnarray}
  \frac {d S}{d t} &=& -\frac{\beta S I}{N}\
  \frac {d E}{d t} &=&\frac{\beta S I}{N}- k E\
   \frac {d I}{d t} &=& k E -\gamma I\
   \frac {d R}{d t} &=& \gamma I \
   \frac {d C}{d t} &=& k E\
\end{eqnarray}
where $C(t)$ is not an epidemiological state but serves to keep track of the cumulative number of Ebola cases from the time of onset of symptoms.
$1/k$ Mean incubation period (days)
$1/\gamma$ Mean infectious period (days)
$\beta=\frac{\beta_0+\beta_1}{2}$
Case for Congo 1995
N=564500
$\beta_0=0.33,\beta_1=0.09$
$1/k=5.3,1/\gamma=5.61$
End of explanation
q=0.5
tI=5
th=tI+np.log(2)/q
beta1+(beta0-beta1)*np.exp(-q*(th-tI))
def bt(b0,b1,tI,q,T,n):
    beta0=b0
    beta1=b1
    dt=T/float(n)
    t=np.linspace(0,T,n+1)
    bI=np.zeros(n+1)
    nI=np.ceil(n*tI/T)
    bI[:nI]=beta0*np.ones(nI)
    bI[nI:]=beta1+(beta0-beta1)*np.exp(-q*(t[nI:]-tI))
    return bI
bI=bt(beta0,beta1,5,0.29,T,n)
plt.xlim(0,100)
plt.plot(t,bI,t,(beta0+beta1)/2.*np.ones(np.size(t)))
def BRN(b):
    s=np.size(b)
    R=np.zeros(s)
    i=0
    for i in np.arange(s):
        J=np.matrix([[-k, b[i]],[k, -gamma]]);
        r=np.max(np.real(la.eigvals(J)));
        R[i]=1+(r*r+(k+gamma)*r)/k/gamma
    return R
R00=BRN(bI)
plt.plot(t,R00)
S,E,I,R=np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)),np.zeros(len(t)) 
S[0]=S0;E[0]=E0;I[0]=I0;R[0]=R0
for i in range(len(t)-1):
    S[i+1] = S[i]-dt*beta*S[i]*I[i]/N
    E[i+1] = E[i]+dt*(beta*S[i]*I[i]/N-k*E[i])
    I[i+1] = I[i]+dt*(k*E[i]-gamma*I[i]) 
    R[i+1] = N-S[i+1]-I[i+1]-E[i+1]
plt.plot(t,S/N,t,I/N,'r--',t,E/N,'y',t,R/N,'k-',t,R00) 
#plt.ylim(0,1)
plt.legend(['Suspective','Infected','Exposed','Recoverd','$R_0$'])
Explanation: Intervention Case
After diease breakout, there ahould be some intervention strategies made to control the spread of Ebola.
The transmission rate, $\beta$, could be modeled reasonably as follows:
$$
\beta ( t ) = \left{\begin{array}{l}
     \beta_{0},\qquad\qquad\qquad\qquad t \leqslant t_{I}\
     \beta_{1} + ( \beta_{0} - \beta_{1} ) e^{-q ( t-t_{I} )},  \quad  t
     \geqslant t_{I}
   \end{array}\right.
$$
i.e. decreases fom $\beta_0$ to $\beta_1$ at $t_I$, the time at which interventions start, and continus to be decreaing exponentially. 
Control of Transition rate, $q$
Another interpretation of the parameter q can be given in terms
 $$t_h=t_I+ln(2)/q$$
 where $t_h$ is the time at which $\beta(t)=\frac{\beta_1+\beta_0}{2}$
End of explanation
from sympy import Symbol , dsolve , Function , Derivative , Eq, lambdify, symbols
x=Symbol('x')
y0=Symbol('y0')
y=Function('y')
f_=Derivative(y(x),x)-y(x)
sol=dsolve(Derivative(y(x),x)-y(x),y(x),ics={y(0):y0})
print(sol)
C1=Symbol('C1')
ysol=sol.subs({C1:y0})
print(ysol)
z= lambdify(x,ysol.rhs)
z= lambdify(x,y0*exp(x))
z(1)
A,N=symbols("A N")
g_=Derivative(y(x),x)-A*y(x)*(N-y(x))
sol2=dsolve(g_,y(x))
sol2
Explanation: Symboilc Solution dor ODE's 
Up to now, we only solve ODE's by numeric schemes (approximately); can we solve them theoretically? Just as well-know fact:
 Not all integrals can be integrated!
Yes, "Definite No"!    
But we still have chance to solve ODE's by "sympy"!
Example (Exponential Growth)
The solution of $y'(t) = y(t), y(0) = y0$ is $y(t)=y0 \exp(t)$
End of explanation
<END_TASK>