1017.sagews

def eulers_method_vector(f,g,x0,y0,h,N,xName,xmin,xmax,yName,ymin,ymax):
    poses = []
    currentPos = (x0,y0)
    for i in xrange(N):
        delta = [X.subs({xName:currentPos[0],yName:currentPos[1]}) for X in (f,g)]
        poses.append(currentPos)
        currentPos = (currentPos[0] + h*delta[0], currentPos[1] + h*delta[1])
        if (currentPos[0] < xmin) or (currentPos[0] > xmax) or (currentPos[1] < ymin) or (currentPos[1] > ymax):
            break
        #if we've looped back, stop because we're misleading
        if (i>50) and (abs(currentPos[0]-poses[0][0]) < 5*h) and (abs(currentPos[1]-poses[0][1]) < 5*h):
            break
    return line2d(poses) 
def phase_portrait(f,g,xName,xmin,xmax,yName,ymin,ymax, points=None):
    lenFG = sqrt(f^2+g^2)
    newF = f/lenFG
    newG = g/lenFG
    VF = plot_vector_field((newF,newG), (xName, xmin, xmax), (yName, ymin, ymax))
    pts = [((xmax-xmin)*random()+xmin, (ymax-ymin)*random()+ymin) for nothing in xrange(10)]
    #pts = [(xmin + ((xmax-xmin)/5)*i, ymin+((ymax-ymin)/5)*j) for i in xrange(5) for j in xrange(5)]
    plots = sum([eulers_method_vector(f,g,X[0],X[1], 0.01, 1000, xName, xmin, xmax, yName, ymin, ymax) for X in pts])
    return VF + plots

var('x y')

(x, y)

phase_portrait((-y-x)/(x^2+y^2), (x-y)/(x^2+y^2), x, -3, 3, y, -3, 3)