import sympy from sympy import * t, A2, A1, A0, B2, B1, B0, C, D, E = symbols('t A2 A1 A0 B2 B1 B0 C D E') Y = ( A2 * t**2 + A1 * t + A0 ) * exp(t) * sin(2*t) + ( B2 * t**2 + B1 * t + B0 ) * exp(t) * cos(2*t) + C * exp(-t) * cos(t) + D * exp(-t) * sin(t) + E * exp(t) pt1 = diff(Y, t, 2) pt2 = 3 * diff(Y, t, 1) pt3 = 2 * Y residual = ( pt1 + pt2 + pt3 - ( exp(t) * ( t**2 + 1 ) * sin(2*t) + 3 * exp(-t) * cos(t) + 4 * exp(t) ) ).expand() print residual # Output is (after regrouping): # # 2*A0*exp(t)*sin(2*t) + 5*A1*exp(t)*sin(2*t) + 2*A2*exp(t)*sin(2*t) - 10*B0*exp(t)*sin(2*t) - 4*B1*exp(t)*sin(2*t) - exp(t)*sin(2*t) # # + 10*A0*exp(t)*cos(2*t) + 4*A1*exp(t)*cos(2*t) + 2*B0*exp(t)*cos(2*t) + 5*B1*exp(t)*cos(2*t) + 2*B2*exp(t)*cos(2*t) # # + 2*A1*t*exp(t)*sin(2*t) + 10*A2*t*exp(t)*sin(2*t) - 10*B1*t*exp(t)*sin(2*t) - 8*B2*t*exp(t)*sin(2*t) # # + 10*A1*t*exp(t)*cos(2*t) + 8*A2*t*exp(t)*cos(2*t) + 2*B1*t*exp(t)*cos(2*t) + 10*B2*t*exp(t)*cos(2*t) # # + 2*A2*t**2*exp(t)*sin(2*t) - 10*B2*t**2*exp(t)*sin(2*t) - t**2*exp(t)*sin(2*t) # # + 10*A2*t**2*exp(t)*cos(2*t) + 2*B2*t**2*exp(t)*cos(2*t) # # - C*exp(-t)*sin(t) - D*exp(-t)*sin(t) # # - C*exp(-t)*cos(t) + D*exp(-t)*cos(t) - 3*exp(-t)*cos(t) # # 6*E*exp(t) - 4*exp(t) #