Why is my MATLAB's bode plot wildly off? I will give out all the details in case it is relevant. I have a MIMO state space system. I find its bode plot using MATLAB and separately using Mathematica. The plot from MATLAB is wildly off compared to the correct bode plot. However, the Mathematica plot is quite close to the correct one. Here's the interesting thing. The transfer function whose bode plot is being taken in both the softwares is calculated to be the same. Have I made some mistake or is it numerical errors contributing to the problem? Here is relevant part of my code:
A=[0,0,0,0,1,0,0,0;
    0,0,0,0,0,1,0,0;
    0,0,0,0,0,0,1,0;
    0,0,0,0,0,0,0,1;
    -297.3,163.5,0,0,0,0,0,0;
    162.9,-267.2,104.2,0,0,0,0,0;
    0,57.8,-74.2,16.4,0,0,0,0;
    0,0,16.4,-16.4,0,0,0,0];
B=[0,0,0,0,0;
    0,0,0,0,0;
    0,0,0,0,0;
    0,0,0,0,0;
    131.4,0.046,0,0,0;
    0,0,0.045,0,0;
    0,0,0,0.025,0;
    0,0,0,0,0.25];
S=[1,0,0,0,0,0,0,0;
    0,0,0,1,0,0,0,0];
D=zeros(2,5);
sys=ss(A,B,S,D)
systf=tf(sys);
s1a1=systf(1,2);
bode(s1a1);

[MATLAB plot which is way off the correct one][1][Mathematica Plot which is quite close to the correct one]2
EDIT: The transfer function being plotted here is:
$$\frac{0.046s^6-5.522\times 10^{-17}s^5+16.46 s^4 - 1.748\times 10^{-14} s^3 + 880.1 s^2 - 7.808\times 10^{-14} s + 7108}{s^8 - 4.373\times 10^{-15} s^7 + 655.1 s^6 - 2.258\times 10^{-12} s^5 + 9.887\times 10^{4} s^4 - 2.144\times 10^{-10} s^3 + 3.43\times 10^{6} s^2 - 4.353\times 10^{-9} s + 2.069\times 10^{7}}$$
I suspect that the reason for the difference is the number of samples used in the plot, since the extreme points are likely going to infinity (or very large numbers) such that if less points are used, those large values are simply missed in the computation. If the OP uses the same number of points in each, then the same result should be achieved- and for all the points shown if they were super-imposed on each plot it would be clear that for the points given, the results are the same.
My first suspicion before the OP added the plots was that it was due to the units of the axis. That is not the case but leaving this part of my answer below as it is another common reason for Bode plot differences:
The units of the frequency axis should be confirmed as they could be in units of:
Continuous domain units of frequency:
Hertz, or cycles/sec
radians/sec
Discrete domain units of (normalized) frequency:
cycles per sample (unique domain is $[-1,1)$)
radians/sample (unique domain is $[-\pi, \pi)$)