#A script to illustrate how to solve systems of ODEs numerically, and to
#draw a direction field for a pair of autonomous ODEs in Sage.

#Consider a pair of ODEs x1' = f(t,x1,x2), x2' = g(t,x1,x2) where
#f(t,x1,x2) = x1-x2^2, g(t,x1,x2) = x1*x2+x1 (these happen to be autonomous,
#but this isn't necessary for the "dsolve_system_rk4()" below.)
def f(t,x1,x2):
    return(x1-x2^2)
def g(t,x1,x2):
    return(x1*x2+x1)

#Numerical Solution: We can solve the system numerically with initial
#data x1(0) = 4, x2(0) = 1 on the interval 0 <= t <= 2 with the commands
var('t x1 x2');
sol = desolve_system_rk4([f(t,x1,x2),g(t,x1,x2)],[x1,x2],ics=[0,4,1],ivar=t,end_points=2,step=0.05)
#The last argument is the stepsize (can be omitted, defaults to 0.1).

#The output "sol" is a list of the form (t,x1,x2) at specified times.
#To isolate the values of x1 and plot, execute
x1vals = [ [i,j] for i,j,k in sol]
plt1 = list_plot(x1vals,plotjoined=True, color='red')
show(plt1)

#Or to also plot x2 and overlay the graph of both x1(t) and x2(t)
x2vals = [ [i,k] for i,j,k in sol]
plt2 = list_plot(x2vals,plotjoined=True, color='blue')
show(plt1+plt2)

#A parametric plot of the solution trajectory can be obtained with
x1x2vals = [ [j,k] for i,j,k in sol]
list_plot(x1x2vals,plotjoined=True, color='blue').show()

#Direction Fields: If the ODEs are autonomous, we can sketch a
#direction field on a range a < x1 < b, c < x2 < d, with or without
#solution trajectories, by using the supplied command "draw_auto_field.sage".
load('draw_auto_field.sage')

#Make sure f and g do not depend on t.
def f(x1,x2):
    return(x1-x2^2)
def g(x1,x2):
    return(x1*x2+x1)

#Construct a direction field on the range -5 <= t <= 5, -5 <= u <= 5.
draw_auto_field(f, g, [-5, 5, -5, 5])

#Construct a direction field with solutions satisfying x1(0) = 4, x2(0) = 1,
#and x1(0) = 3, x2(0) = -2. Use stepsize 0.01 (last argument) to sketch
#trajectories.
draw_auto_field(f, g, [-5, 5, -5, 5], [[4,1],[3,-2]],0.01)