How to compute Laplace Transform in Python?
Update : Given the approach started in the OP's Github code I have this suggestion :
Observe that the unilateral Laplace Transform given as:
$$X(s) = \int_0^\infty x(t)e^{-st}dt$$
Is just the Fourier Transform of a causal function with a weighting exponential:
$$X(s) = \int_0^\infty x(t)e^{-(\sigma+j\omega)t}dt$$
$$X(s) = \int_0^\infty [e^{-\sigma t}x(t)]e^{-j\omega t}dt$$
Which is the Fourier Transform of $e^{-\sigma t}x(t)$.
So to proceed, learn how to produce surface plots in python, and then index through $\sigma$ within the Region of Convergence (see below) and compute the FFT of $e^{-\sigma t}x(t)$  to create your complex surface values given $\sigma$ and $\omega$.

The challenge with what you are trying to do is that the Laplace Transform is a function of the complex variable &quot;s&quot;, so for each possible value of &quot;s&quot; (which is simply the set of all complex numbers) the Laplace Transform would have a complex result with a magnitude and phase.
So when you try to plot this, you could plot the magnitude and phase as two separate plots but these would be 3D surface plots showing either the result of the magnitude or the phase over a 2D surface, similar otherwise to how we plot the magnitude and phase of the Fourier Transform as a 2D function over the line representing all frequencies.
For the graphical representation of the Laplace Transform, we typically just show the locations where that function goes to infinity (poles) or is zero (zeros). In fact every other sample on the surface is uniquely defined by the pole and zero locations, so that is all we need to show to define it.
This is best shown with an example. Since I already have the graphic for this particular case, consider the time domain function of a decaying sinusoid given by the formula below and the plot below that where we see in the dashed red line the envelope for the decaying function $2e^{-0.2t}$. $u(t)$ is the step function which is $0$ for time $t&lt;0$ and $1$ for time $t \ge 0$.
$$ x(t) = u(t)2e^{-0.2t}sin(0.5t)$$

To get the Laplace Transform (easily), we decompose the function above into exponential form and then use the fundamental transform for an exponential given as :
$$\mathscr{L}\{u(t) e^{-\alpha t}\} = \frac{1}{s+\alpha}$$
This is the unilateral Laplace Transform (defined for $t = 0$ to $\infty$), and this relationship goes a long way since we can describe the response of any causal linear system using such exponential forms.
So the equation above, assuming $t&gt;0$, and using Euler's identity becomes :
$$ x(t) = u(t)2e^{-0.2t}sin(0.5t)$$
$$ = 2e^{-0.2t}\frac{(e^{+j0.5t}-e^{-j0.5t})}{2j}$$
$$ = je^{-0.2t}(e^{+j0.5t}-e^{-j0.5t})$$
$$ = je^{-(0.2-j0.5)t} -je^{-(0.2+j0.5)t}$$
Which we can then easily take the Laplace Transform to get a function of the complex variable $s$ :
$$ X(s) = \frac{j}{s+0.2+j0.5}-\frac{j}{s+0.2-j0.5} = \frac{1}{(s+0.2+j0.5)(s+0.2-j0.5)}$$
The graph of the magnitude of this is shown below, which is a surface plot for all values of $s$, where $s$ is a complex variable with real and imaginary components traditionally described as $s = \sigma +j\omega$.  So we simply choose a particular complex value $s$, and plot the magnitude of the result of $X(s)$.

Caveat : This is indeed the magnitude of the Laplace Transform if the Laplace Transform exists. For causal functions, the Laplace Transform does not exist for all values of s whose real part is to the left of the leftmost pole. Why? Because for any values of $s$ in that region (left of the leftmost pole), the Laplace Transform (given by the integral) will not converge (grows to infinity). While for all values of s at and to the right of the leftmost pole, the transform will converge to be the magnitude on the surface plotted above. So the proper graphic would only be valid to the right of the right most pole. This can be confusing as $X(s)$ certainly converges as was done to make this plot, but it is in the transformation itself that the result cannot be obtained.
This is very clear if we consider the envelope in our example function $x(t)$ which was given by $2e^{-0.2t}$ for all $t&gt;0$, and the Laplace Transform $X(s)$ for $s = -1$ :
$$X(s) = \int_0^\infty x(t)e^{-st}dt$$
$$X(s= -1) = \int_0^\infty 2e^{-0.2t}e^{t}dt = \int_0^\infty 2e^{0.8t}dt$$
which will not converge since the function $e^{0.8t}$ continuously grows larger for larger $t$.
Hopefully after reading this you no longer feel the need to plot the Laplace Transform, and in practical application a plot of it is never used beyond showing the pole and zero locations. However what is very useful is knowing that the Fourier Transform is the Laplace Transform when $s = j\omega$. You can see this if you compare the two equations, and the small breakout in upper right-hand corner of the plot above is also showing this, which is the Frequency Response specifically. The one precaution is that the Fourier Transform is often given as a bilateral function (t extending from $-\infty$ to $\infty$) so to be truly equivalent unless the function is declared to be causal, we must be using the bilateral Laplace Transform for the two to be exactly identical (which is also seldom used). This explains why the Fourier Transform of a sine wave would appear as two impulses at $\pm \omega_c$ while as you might deduce from the plot above, the Laplace Transform of the sine wave would appear more as tent poles on the $j\omega$ axis. You would see this as well with the Fourier Transform if you did the transform of $u(t)cos(\omega t)$ specifically.