# Revision history [back]

Here is one solution using ode_solver.

var('S,E,I,W')

S0=2600
E0=1
I0=1
W0=200
bR= 37
k= 8760
m= .0142
r= .9
B= .0255*(1/W0)
A= 1
Y=100000
mW=26
timedata=srange(0,10,0.05)

f=[bR-B*S*W-m*S+k*I,B*S*W-A*E-m*E,A*E-k*I-m*I,Y*I-mW*W]
soln=desolve_odeint(f,[S0,E0,I0,W0],timedata,[S,E,I,W])

splot=line(zip(timedata,[s[0] for s in soln]),color='green')
eplot=line(zip(timedata,[s[1] for s in soln]),color='blue')
iplot=line(zip(timedata,[s[2] for s in soln]),color='red')
wplot=line(zip(timedata,[s[3] for s in soln]),color='black')
show(splot)
show(eplot)
show(iplot
show(wplot)


Here is one solution using ode_solverdesolve_odeint.

var('S,E,I,W')

S0=2600
E0=1
I0=1
W0=200
bR= 37
k= 8760
m= .0142
r= .9
B= .0255*(1/W0)
A= 1
Y=100000
mW=26
timedata=srange(0,10,0.05)

f=[bR-B*S*W-m*S+k*I,B*S*W-A*E-m*E,A*E-k*I-m*I,Y*I-mW*W]
soln=desolve_odeint(f,[S0,E0,I0,W0],timedata,[S,E,I,W])

splot=line(zip(timedata,[s[0] for s in soln]),color='green')
eplot=line(zip(timedata,[s[1] for s in soln]),color='blue')
iplot=line(zip(timedata,[s[2] for s in soln]),color='red')
wplot=line(zip(timedata,[s[3] for s in soln]),color='black')
show(splot)
show(eplot)
show(iplot
show(wplot)