PyStarshade simulates high-contrast direct imaging of exoplanets with starshades. It models the optical performance of starshade-based missions by propagating complex electric fields through three planes (starshade, telescope aperture, and focal plane) using Fresnel or Fraunhofer diffraction computed with Bluestein FFTs. The toolkit computes diffracted fields and point-spread functions for a variety of starshade and telescope configurations to evaluate metrics such as core throughput, contrast, and inner working angle, and can simulate imaging for discretized exoplanet scenes while interfacing with HCIPy (ascl:2602.012) to generate telescope apertures.