237.sagews

#
# dy/dx = 0.1*sqrt(y) + 0.4*x**2
#

from pylab import *

def fn(x, y):
	return 0.1*sqrt(y) + 0.4*x**2

def isocline_fn(c, x):
	return ((c - 0.1*x**2)/(0.4))**2

def xcomp(fn):
	return 1 / sqrt(1 + fn**2)

def ycomp(fn):
	return fn / sqrt(1 + fn**2)

xr = arange(-2,3,0.3)
yr = arange(0,5,0.3)

figure()
X = zeros([len(xr), len(yr)])
Y = zeros([len(xr), len(yr)])
U = zeros([len(xr), len(yr)])
V = zeros([len(xr), len(yr)])
for x in range(0, len(xr)):
	for y in range(0, len(yr)):
                X[x][y] = xr[x]
                Y[x][y] = yr[y]
		U[x][y] = xcomp(fn(xr[x], yr[y]))
		V[x][y] = ycomp(fn(xr[x], yr[y]))

Q = quiver(X, Y, U, V)

# isocline
xr = arange(-2,3,0.1)
for c in arange(0,10,0.05):
	yy = isocline_fn(c, xr)
	plot(xr, yy)

# exact sln

axis([-2,3,0,5])

savefig("df-1.png")