#Parameters
a = .02
b = .05
c = -.03 
d = -.01
# Initial condition
RInit = 1
JInit = 1


# End time
tMax = 20


# Standard SIR model
def ODE_RHS(t, Y):
    (R,J) = Y
    
    dR = a*R + b*J
    dJ = c*R + d*J
   

    return (dR, dJ)

# Set up numerical solution of ODE
solver = ode_solver(function = ODE_RHS,
                    y_0 = (RInit, JInit),
                    t_span = (0, tMax),
                    algorithm = 'rk8pd')

# Numerically solve
solver.ode_solve(num_points = 500)

# Plot solution
show(
     plot(solver.interpolate_solution(i = 0), 0, tMax, legend_label = 'R(t)', color = 'green')
   + plot(solver.interpolate_solution(i = 1), 0, tMax, legend_label = 'J(t)', color = 'red')
) 

# code from https://sage.math.clemson.edu:34567/home/pub/161/
# Thanks to Jan Medlock