def draw_auto_field(f, g, viewrange, *args):
    #Usage: draw_auto_field(f, g, viewrange, *args)
    #Draws direction field for x1' = f(t,x1,x2), x2' = g(t,x1,x2)
    #on the range a < x < b, c < y < d where viewrange=[a,b,c,d].
    #If args is defined then ic = args[0]
    #and ic is an optional list of initial conditions for solution curves,
    #at time t = 0, of the form ic = [[a0,b0],[a1,b1],...,[an,bn]]

    var('t x1 x2')

    #Pick off low, high in each variable, make sure all are floating point.
    x1low = viewrange[0].n(); x1high = viewrange[1].n();
    x2low = viewrange[2].n(); x2high = viewrange[3].n();

    #Plot vector field itself with built-in Sage command
    p = plot_vector_field((f(x1,x2),g(x1,x2)), (x1,x1low,x1high), (x2,x2low,x2high))

    #If initial points are supplied, compute solution trajectory through those points.
    if len(args)>0:
        ic = args[0]
        N = len(ic)
        R = 2
        h = args[1]

        for i in range(0,N):
            x1_0 = ic[i][0].n();
            x2_0 = ic[i][1].n();

            sol = desolve_system_rk4([f(x1,x2),g(x1,x2)],[x1,x2],ics=[0,x1_0,x2_0],ivar=t,end_points=R,step=h);
            x1x2vals = [ [j,k] for i,j,k in sol]
            p2 = list_plot(x1x2vals,plotjoined=True, color='red')
            p = p + p2;
            sol = desolve_system_rk4([f(x1,x2),g(x1,x2)],[x1,x2],ics=[0,x1_0,x2_0],ivar=t,end_points=-R,step=h);
            x1x2vals = [ [j,k] for i,j,k in sol]
            p2 = list_plot(x1x2vals,plotjoined=True, color='red')
            p = p + p2;

    show(p,xmin=x1low,xmax=x1high,ymin=x2low,ymax=x2high)

    return