Slightly more Sage-ified version of the above very nice solution:

import scipy.integrate
a=1.0
b=2.0

def fun(t):
    if t<=-b:
        return -a
    elif f<b:
        return t*a/b
    else:
        return a

g=lambda t:fun(t)

N=100
time_step=0.1
time_end=10.0
t0=0.0
x0=[[0.5*k,0.5*k] for k in range(-10,10)]

def f(x,t):
    return [x[1],-x[0]+g(x[0])]

time_range=[t0..time_end, step=time_step]
x0=[[0.5*k,0.5*k] for k in range(-10,10)]
sol_lines = Graphics()
for n in range(10):
    sol = scipy.integrate.odeint(f,x0[n],time_range)
    sol_lines += line(sol,rgbcolor=hue(.3+n/15.0))

x0,x1=var('x0 x1')
p=plot_vector_field ((x1,-x0+g(x0)),(x0,-9,9),(x1,-7,7))

show(sol_lines + p, figsize = [9,7])

-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org

Reply via email to