Differential Plotter

import numpy as np
import matplotlib.pyplot as plt

# Tells the program how big the window should be, and how far apart
# the derivatives should render
xDets={
    "max":10,
    "min":-10,
    "step":1
    }
yDets={
    "max":210,
    "min":-10,
    "step":10
    }
maxSteps=30000
solStepSize=0.01
H=120
initialPoints=[[0,42],[0,101],[0,150]]
colors=['g','m','y','r']
#the derivative function used
def diff(x,y):
    return y*(2-0.01*y)-H

#computes how many times the step goes into the specified range, and rounds up
#to be as close as possible
def getSize(xDets,yDets):
    xSteps=int(np.ceil((xDets["max"]-xDets["min"])/xDets["step"])+1)
    ySteps=int(np.ceil((yDets["max"]-yDets["min"])/yDets["step"])+1)
    return xSteps,ySteps

xSteps,ySteps=getSize(xDets,yDets)
# Establishes arrays for the x and y axis with each entry at the gridlines
x = np.linspace(xDets["min"],xDets["max"],xSteps)
y = np.linspace(yDets["min"],yDets["max"],ySteps)

# creates a pair for every point so X, Y can represent the whole grid
X, Y= np.meshgrid(x,y)

# Computes the derivative at all points
Diffs=diff(X,Y)

# Constructs a normalized vector in the direction of the derivative
Lens=np.sqrt(Diffs**2+1)
U=1/Lens
V=Diffs/Lens


#plots it
fig, ax=plt.subplots()
fig.savefig('graph.png')
ax.quiver(X,Y,U,V,color="C0",angles="xy",
          scale_units="xy",width=0.0015)

def plotSoln(x0,y0):
    points = []
    x,y=x0,y0
    steps=0
    while steps<maxSteps and isInBounds(x,y)and [x,y] not in points:
        diffs=[1,diff(x,y)]
        points.append([x,y])
        steps+=1
        norm=np.sqrt(diffs[0]**2+diffs[1]**2)
        x+= solStepSize*diffs[0]/norm
        y+= solStepSize*diffs[1]/norm
    steps=0
    x,y=x0,y0
    diffs =[1,diff(x,y)]
    steps+=1
    norm=np.sqrt(diffs[0]**2+diffs[1]**2)
    x-= solStepSize*diffs[0]/norm
    y-= solStepSize*diffs[1]/norm
    while steps<maxSteps and isInBounds(x,y)and [x,y] not in points:
        points.insert(0,[x,y])
        diffs =[1,diff(x,y)]
        steps+=1
        norm=np.sqrt(diffs[0]**2+diffs[1]**2)
        x-= solStepSize*diffs[0]/norm
        y-= solStepSize*diffs[1]/norm
    return points

def isInBounds(x,y):
    return xDets["min"]<=x and xDets["max"]>=x and yDets["min"]<=y and yDets["max"]>=y
for idx,pt in enumerate(initialPoints):
    soln=np.array(plotSoln(pt[0],pt[1]))
    plt.plot(soln[:,0],soln[:,1],color=colors[idx])

ax.set(xlim=(xDets["min"],xDets["max"]),ylim=(yDets["min"],yDets["max"]))
ax.axhline(y=0, color='k')
ax.axvline(x=0, color='k')

plt.show()