Revision history [back]

Here's an option using ode_solve.

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

T = ode_solver()
T.function=lambda t,y:[bR-B*y[0]*y[3]-m*y[0]+k*y[2],B*y[0]*y[3]-A*y[1]-m*y[1],A*y[1]-k*y[2]-m*y[2],Y*y[2]-mW*y[3]]
T.ode_solve(y_0=[S0,E0,I0,W0],t_span=[0,10],num_points=500)


a=T.solution
tssoln=[[soln[0],soln[1][0]] for soln in a]
tesoln=[[soln[0],soln[1][1]] for soln in a]
tisoln=[[soln[0],soln[1][2]] for soln in a]
twsoln=[[soln[0],soln[1][3]] for soln in a]

show(line(tssoln,color='green'))
show(line(tesoln,color='blue'))
show(line(tisoln,color='red'))
show(line(twsoln,color='black'))