This is the same question as this one, except for Python instead of Mathematica. Basically, the MATLAB software PPLANE is a staple in ODE courses. Is there a Python equivalent?

I don't know much about the software since my class just started using it, but I basically need it for plotting nullclines and trajectories, as well as finding fixed points and their stability. PPLANE can basically do all of this through a nice GUI without having to write any code.

Sample outputs from pplane:

PPLANE example slope field

PPLANE example slope field

  • 1
    $\begingroup$ What does the software do? $\endgroup$
    – nicoguaro
    Commented Oct 18, 2021 at 3:36
  • 2
    $\begingroup$ following up on @nicoguaro's question, what parts of the software do you need? Plotting some streamlines and fixed points is relatively easy to do yourself. Finding closed orbits is a bit harder. $\endgroup$ Commented Oct 18, 2021 at 6:23
  • $\begingroup$ @ThijsSteel see edit to the question $\endgroup$ Commented Oct 18, 2021 at 14:05
  • $\begingroup$ @nicoguaro see edit ^ $\endgroup$ Commented Oct 18, 2021 at 14:06
  • $\begingroup$ Please point your teaching assistant to the following resource: github.com/MathWorks-Teaching-Resources/…. The old pplane is now severely outdated and no longer works for the latest versions of Matlab. $\endgroup$ Commented Oct 18, 2021 at 14:52

3 Answers 3


I'm not aware of any alternatives written in python.

However, i don't think an alternative is strictly necessary. Firstly, even if you don't know how to write Matlab code, you should be able to figure out the GUI. Secondly, The app https://github.com/MathWorks-Teaching-Resources/Phase-Plane-and-Slope-Field is the latest version of pplane and is written as a Matlab app, so it should be possible to run it without Matlab installed (if you can't get it working without a license, you can maybe ask your TA to help you create an executable).

If neither of these work, you'll probably have to write some python code to simulate the system and plot it yourself.

  • $\begingroup$ There are multiple reasons to want an alternative. You could want to enable your students to run it on their own laptop and use It in homework for example. In this case, a commercial Matlab license is not an option. $\endgroup$
    – BlaB
    Commented Oct 19, 2021 at 10:24
  • $\begingroup$ I believe that it should be possible to run that GUI without a matlab license. But sure, a python/c++ version would be useful to not have to go through the process. You might also want students to be able to edit the graph easily without knowing matlab. $\endgroup$ Commented Oct 19, 2021 at 11:32

I do not really know how much this can help you, but maybe you can use this code for rough sketches of the dynamics of planar vector fields. I know that it does not have the functionality that you may really wish for, but it may come in handy to guide the analysis in the right direction.

import numpy as np
import matplotlib.pyplot as plt

def plot_dynamics(vector_field, x_left, x_right, x_res, y_down, y_up, y_res):
    x, y = np.meshgrid(np.linspace(x_left, x_right, x_res), np.linspace(y_down, y_up, y_res))
    Vx, Vy = vector_field(x, y)
    if type(Vx) != object:
        Vx = Vx * np.ones(x.shape, dtype=float)
    if type(Vy) != object:
        Vy = Vy * np.ones(x.shape, dtype=float)
    fig, ax = plt.subplots()
    #ax.set_aspect( 1 )
    ax.streamplot(x, y, Vx, Vy)
    return None

# type the formulas for the x and y components of the vector fields 
# (use np.cos and np.sin etc if not polynomial vector fields):
def V(x, y):    
    return  ( 2*x - y + 3*(x**2-y**2) + 2*x*y,   x - 3*y - 3*(x**2-y**2) + 3*x*y )

def f(x, y):    
    return  ( 1,  y**2 - x ) 

plot_dynamics(V, -2, 4, 100, -4, 2, 100)
plot_dynamics(f, -2, 10, 100, -4, 4, 100)

enter image description here

enter image description here


I would like to provide with my codes with a result resembling that of matlab, although not 100% the same. Codes are maintained and extended to more general cases via pplane for python.

# -*- coding: utf-8 -*-
pplane for phase portrait
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator
from scipy.optimize import fsolve

def eqnXDotYDot(x, y): 
    dx = 2*x -y +3*(x**2-y**2) +2*x*y
    dy = x-3*y-3*(x**2-y**2)+3*x*y
    return  ( dx, dy )

def func(variables):
    x, y = variables
    dx, dy = eqnXDotYDot(x, y)
    return (dx, dy)

def fixedPoints(X, Y):
    tolerance = 1e-8
    fPtX = []
    fPtY = []
    for i in range(nx):
        for j in range(ny):
            xFixedPoint, yFixedPoint = fsolve(func,(X[i,j],Y[i,j]))
            if all(np.isclose(eqnXDotYDot(xFixedPoint, yFixedPoint), [0.0, 0.0])) == True:
                if i == 0 and j == 0:
                    oldX, oldY = xFixedPoint, yFixedPoint
                    if abs(oldX - xFixedPoint)<tolerance and abs(oldY - yFixedPoint)<tolerance:
                    oldX, oldY = xFixedPoint, yFixedPoint
    return (fPtX,fPtY)

# Grid of x, y points
nx, ny = 100, 100
minX, maxX = -2, 4
minY, maxY = -4, 2

x = np.linspace(minX, maxX, nx)
y = np.linspace(minY, maxY, ny)
X, Y = np.meshgrid(x, y)    

# field vector
dx, dy = eqnXDotYDot(X,Y)
# plot phase portrait with vector field
plt.title("Phase Portrait")
ax = plt.gca()
#fig, ax = plt.subplots()
minorLocatorX = AutoMinorLocator(2) # number of minor intervals per major inteval
minorLocatorY = AutoMinorLocator(2)
ax.xaxis.set_minor_locator(minorLocatorX) # add minor ticks on x axis
ax.yaxis.set_minor_locator(minorLocatorY) # add minor ticks on y axis

speed = np.sqrt(dx**2 + dy**2)
lw = speed  / speed.max()

ax.streamplot(x, y, dx, dy, linewidth=lw, density=1,color='b', arrowstyle='-',broken_streamlines=False)

ax.contour(X,Y,dx,levels=[0], linewidths=1, colors='r')
ax.contour(X,Y,dy,levels=[0], linewidths=1, colors='y')

fPtX,fPtY = fixedPoints(X, Y)
ax.scatter(fPtX,fPtY,s = 80, facecolors='none',edgecolors='g')

nx, ny = 15, 15
x = np.linspace(minX, maxX, nx)
y = np.linspace(minY, maxY, ny)
X, Y = np.meshgrid(x, y)  
dx, dy = eqnXDotYDot(X,Y) 
q = ax.quiver(X, Y, dx, dy,color='g',angles = 'uv', headlength=2,headaxislength=2)

ax.set_xlim(1.05*minX, 1.05*maxX)
ax.set_ylim(1.05*minY, 1.05*maxY)

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.