Why so complicated?
Why can’t the code be as simple as:
# Random code I wrote just to illustrate the issue - won't work
from sympy import symbols, solve
from sympy.integrals.transforms import laplace_transform, inverse_laplace_transform
s = symbols('s')
t = symbols('t', positive = True)
y = Function('y')(t)
L = Function('L')(s)
# input
y0 = [2, 1]
deq = y.diff(t,2) + y.diff(t,1) - 2*y - 4
algeq = laplace_transform(eq, s,t, noconds = True)
algsoln = solve(algeq, L)[0]
soln = inverse_laplace_transform(algsoln, s, t, noconds = True)
print("Differential Equation: ", deq, "=0")
print("Solution: y(t) = ", soln)
The trouble arises when you try to find the Laplace Transform of derivative of y, which laplace_trasform function cannot handle yet. (Try it yourself – print(laplace_transform(y.diff(t), t, s))
So we had to take a shortcut. Since it was supposed to be taught to students and it had to be as simple as possible (we are allowed to cheat a bit on code’s integrity), we settled for the following, which inputs the transformed equation instead of the differential equation itself. The friend who proposed it came up with this code (with minor touches from me for better readability):
from sympy import *
from sympy.integrals.transforms import laplace_transform, inverse_laplace_transform
s, L = symbols('s, L')
t= symbols ('t', positive =True)
y0 = 2
y10 = 1 # storing y(0) and y’(0)
Ly2 = s**2*L-s*y0-y10
Ly1= s*L-y0
Ly = L
algeq = Eq(Ly2 + Ly1 - 2*Ly, laplace_transform(4, t, s, noconds = True))
algsoln = solve(algeq, L)[0]
soln = inverse_laplace_transform(algsoln, s, t, noconds = True)
print("Solution: y(t) = ", soln)
We decided to settle for this until we came to the third problem: 
Mathematicos 