You are approximating a definite integral with cumtrapz it won't give you the same result as the integrated equation unless you add a constant and plot with the given x coordinates:
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as it
x = np.arange(-10,10, 0.01) # start,stop,step
f = x**2
So why don't you very literally do what is being asked.
The difference between two curves is the surface area between the two curves. This is 1 number, which you can then plot as a function of the number of terms in the Taylor expansion. There is however a small snag, the area will depend on the range of x's you consider. So there are two options, either you ...
Sounds like you want to omit the use of a MATLAB toolkit/Octave package. Of course you can build the bode diagram from scratch:
% transfert function as anonymous function
h = @(s)(1)./(2*s.^2+3*s+4);
% frequency vector
w = 2*pi*logspace(-2,2,1000);
% magnitude & phase estimation
mag = abs(h(1j*w));
magDB = 20*log10(mag);
phaseDeg = rad2deg(angle(h(1j*w))...