The article presents programs for the approximate calculation of Lebesgue–Stieltjes integrals in Python, which are not currently available in the SymPy and SciPy packages. Those packages include only functions for calculating single and multiple Riemann integrals. To implement these programs, there has been made the correction of classical mathematical expressions, which determine the Lebesgue–Stieltjes integrals, and synthesized algorithms suitable for the development of programs for the approximate calculation of these integrals in Python. The feature of the algorithm synthesized for the approximate calculation of the Lebesgue integral is that the Lebesgue measure of a discrete function given on a zero-dimensional set of points located on the segment of its argument is a monotonic continuous function of the coordinate of the functional axis. This axis value increases from zero at the point of the minimum value of this function to a value equals to the length of the segment of the functional axis in the range from the minimum value of this function to its maximum value. In this algorithm, the values of a discrete Lebesgue-integrated function are sorted to form an ascending sequence, the measure of each value of which is given by a segment of the functional axis within adjacent values of this sequence in the direction of its growth. The developed Python programs for Lebesgue–Stieltjes integration contain standard already known program functions of this programming language. The article shows that the proposed programs can be useful for scientists who work on problems of systems analysis with discrete models.