We present PyBDR, a Python reachability analysis toolkit based on set-boundary analysis, which centralizes on widely-adopted set propagation techniques for formal verification, controller synthesis, state estimation, etc. It employs boundary analysis of initial sets to mitigate the wrapping effect during computations, thus improving the performance of reachability analysis algorithms without significantly increasing computational costs. Beyond offering various set representations such as polytopes and zonotopes, our toolkit particularly excels in interval arithmetic by extending operations to the tensor level, enabling efficient parallel interval arithmetic computation and unifying vector and matrix intervals into a single framework. Furthermore, it features symbolic computation of derivatives of arbitrary order and evaluates them as real or interval-valued functions, which is essential for approximating behaviours of nonlinear systems at specific time instants. Its modular architecture design offers a series of building blocks that facilitate the prototype development of reachability analysis algorithms. Comparative studies showcase its strengths in handling verification tasks with large initial sets or long time horizons. The toolkit is available at https://github.com/ASAG-ISCAS/PyBDR .