from sympy import *
init_printing()
omega = symbols('omega')
A, sigma, omega_0 = symbols('A sigma omega_0', real = True, positive = True)

C = symbols('C', real = True, positive = True)
I_m = Integral(A*exp(-C*((omega-omega_0)/sigma)**2), (omega, 0, oo))
I_p = Integral(A*exp(-C*((omega+omega_0)/sigma)**2), (omega, 0, oo))
I_m-I_p

x = symbols('x')
I_m = I_m.transform((omega-omega_0)/sigma, x)
I_p = I_p.transform((omega+omega_0)/sigma, x)
I_m-I_p

2*A*sigma*Integral(exp(-C*x**2), (x, 0, omega_0/sigma))

_.doit()

epsilon_0 = symbols('epsilon_0', real = True, positive = True)
NumCoeff = A*epsilon_0*sigma
SW_m = Integral(NumCoeff*exp(-C*((omega-omega_0)/sigma)**2), (omega, 0, oo))
SW_p = Integral(NumCoeff*exp(-C*((omega+omega_0)/sigma)**2), (omega, 0, oo))
SW_m-SW_p

SW_m = SW_m.transform((omega-omega_0)/sigma, x)
SW_p = SW_p.transform((omega+omega_0)/sigma, x)
SW_m-SW_p

SW_1 = Integral(omega_0*exp(-C*x**2), (x, -omega_0/sigma, oo))
SW_2 = Integral(omega_0*exp(-C*x**2), (x, omega_0/sigma, oo))
SW = SW_1 + SW_2
SW

NumCoeff *= omega_0
SW_1 = Integral(exp(-C*x**2), (x, -omega_0/sigma, omega_0/sigma))
SW_2 = Integral(exp(-C*x**2), (x, omega_0/sigma, oo))
SW = SW_1 + 2*SW_2
SW

SW = SW.subs(C, 4*ln(2))
result = SW.doit().simplify()
result *= NumCoeff 
result