Nonlinear Optical Spectroscopy is a well-developed field with theoretical and experimental advances that have benefited multiple disciplines, including chemistry, biology, and physics. However, for the accurate interpretation of the corresponding multi-dimensional spectra, there is a need for precise quantum dynamical simulations based on model Hamiltonians. In this article, we present the initial release of our code, QuDPy (Quantum Dynamics in Python), which provides a robust numerical platform for performing quantum dynamics simulations based on model systems, including open quantum systems. A distinguishing feature of our approach is the ability to specify various high-order optical response pathways in the form of double-sided Feynman diagrams through a straightforward input syntax. This syntax outlines the time-ordering of ket-sided or bra-sided optical interactions acting on the time-evolving density matrix of the system. We utilize the quantum dynamics capabilities of QuTip to simulate the spectral response of complex systems, allowing us to compute virtually any n-th order optical response of the model system. To illustrate the utility of our approach, we provide a series of example calculations. Program summaryProgram Title:QuDPyCPC Library link to program files:https://doi.org/10.17632/5xm9pm24cz.1Developer's repository link:https://github.com/sa-shah/QuDPyLicensing provisions: MIT LicenseProgramming language: Python (v3.7)Supplementary material: Available as Google Colab Files.Example 1: https://tinyurl.com/y3j5jmmrExample 2: https://tinyurl.com/37vwntn5External packages:•QuTip (v.4.7) and dependencies i.e. Numpy, Matplotlib (https://qutip.org/)•UFSS Automatic Diagram Generator (https://github.com/peterarose/ufss)Nature of problem: Accurate quantum simulations of complex systems are required in order to understand and interpret multi-dimensional ultrafast spectroscopic signals. This code provides an open-source/multi-platform method that facilitates the generation of higher-order non-linear optical responses for an arbitrary molecular or material system given a model input Hamiltonian and bath model.Solution method: We use the double-sided Feynman diagram method [1, 2] to generate (symbolically) a set of response functions corresponding to the nth order non-linear response of the system to a series of laser pulses using the UFSS package [3] We then perform a series of accurate quantum dynamics calculations using the QuTip package [4] to generate the numerical response and spectra which correspond to specific experimental conditions.