How to plot 3-D bode from the derived transfer function?

I want to plot my derived transfer function in this format as shown in figure below. The transfer function is derived using state matrices and using this following expression: G_vd=q*(inv((s*I)-A))*f

I used symb to get the transfer function with 'ws=2*pi*Fs - switching frequency' as a variable. 'wo=2*pi*Fo' is a constant and I can vary 'ws' as per the ratio needed. So, now I have a transfer function with variable 'ws' and the complex frequency 's'.

This is what I obtained after simplifying. How to proceed further? I found answers where manual multiplication was performed. Also, I want to use surf command to obtain a neat surface.

G_vd =-(1.6872e+30*ws*(2.0616e+11*ws^2 - 1.5981e+19)*(- 1.0434e+209*s^14*ws^2 + 9.7644e+213*s^14*ws + 3.4484e+220*s^14 - 1.2991e+213*s^13*ws^2 + 1.2158e+218*s^13*ws + 4.2937e+224*s^13 - 5.2168e+209*s^12*ws^4 + 3.9142e+214*s^12*ws^3 + 6.8734e+220*s^12*ws^2 + 1.2884e+226*s^12*ws + 3.4203e+232*s^12 - 5.6241e+213*s^11*ws^4 + 4.8516e+218*s^11*ws^3 + 9.9349e+224*s^11*ws^2 + 9.4516e+229*s^11*ws + 2.8572e+236*s^11 - 1.0434e+210*s^10*ws^6 + 5.8925e+214*s^10*ws^5 + 3.3573e+220*s^10*ws^4 + 4.185e+226*s^10*ws^3 + 6.8393e+232*s^10*ws^2 + 3.2035e+237*s^10*ws + 1.131e+244*s^10 - 9.5051e+213*s^9*ws^6 + 7.2633e+218*s^9*ws^5 + 5.6111e+224*s^9*ws^4 + 3.0823e+230*s^9*ws^3 + 7.0634e+236*s^9*ws^2 + 9.3038e+240*s^9*ws + 4.8152e+247*s^9 - 1.0434e+210*s^8*ws^8 + 3.9566e+214*s^8*ws^7 - 1.3383e+219*s^8*ws^6 + 4.8276e+226*s^8*ws^5 + 3.4139e+232*s^8*ws^4 + 9.5864e+237*s^8*ws^3 + 2.2626e+244*s^8*ws^2 + 3.8015e+246*s^8*ws + 1.2478e+255*s^8 - 7.7621e+213*s^7*ws^8 + 4.8387e+218*s^7*ws^7 - 5.7637e+222*s^7*ws^6 + 3.5915e+230*s^7*ws^5 + 4.2122e+236*s^7*ws^4 + 2.7369e+241*s^7*ws^3 + 1.4076e+248*s^7*ws^2 + 2.6339e+249*s^7*ws + 2.0572e+257*s^7 - 5.2168e+209*s^6*ws^10 + 1.0103e+214*s^6*ws^9 - 3.5371e+220*s^6*ws^8 + 2.2567e+226*s^6*ws^7 + 3.4001e+232*s^6*ws^6 + 9.5471e+237*s^6*ws^5 + 6.907e+241*s^6*ws^4 + 9.913e+246*s^6*ws^3 + 3.7354e+255*s^6*ws^2 + 8.413e+254*s^6*ws + 2.891e+263*s^6 - 3.0095e+213*s^5*ws^10 + 1.2146e+218*s^5*ws^9 + 1.3223e+224*s^5*ws^8 + 1.7168e+230*s^5*ws^7 - 1.3395e+236*s^5*ws^6 + 2.6851e+241*s^5*ws^5 + 1.3718e+248*s^5*ws^4 + 2.9996e+249*s^5*ws^3 + 3.9942e+257*s^5*ws^2 + 2.5473e+257*s^5*ws + 4.5819e+265*s^5 - 1.0434e+209*s^4*ws^12 + 8.4668e+211*s^4*ws^11 - 6.9202e+220*s^4*ws^10 + 3.2845e+225*s^4*ws^9 + 6.818e+232*s^4*ws^8 + 3.146e+237*s^4*ws^7 - 2.2528e+244*s^4*ws^6 + 9.4169e+246*s^4*ws^5 + 3.7335e+255*s^4*ws^4 + 8.9376e+254*s^4*ws^3 + 2.891e+263*s^4*ws^2 + 4.474e+262*s^4*ws + 1.7485e+271*s^4 - 4.276e+212*s^3*ws^12 + 3.47e+215*s^3*ws^11 + 1.3494e+224*s^3*ws^10 + 2.6238e+229*s^3*ws^9 - 1.3454e+236*s^3*ws^8 + 8.7872e+240*s^3*ws^7 + 4.4564e+247*s^3*ws^6 + 1.2063e+249*s^3*ws^5 + 1.9823e+257*s^3*ws^4 + 2.6028e+256*s^3*ws^3 + 9.8906e+262*s^3*ws^2 + 8.3904e+264*s^3*ws + 2.7391e+273*s^3 - 3.4492e+220*s^2*ws^12 + 2.799e+223*s^2*ws^11 + 3.4125e+232*s^2*ws^10 - 1.8193e+235*s^2*ws^9 - 1.1279e+244*s^2*ws^8 + 3.3039e+246*s^2*ws^7 + 1.245e+255*s^2*ws^6 + 2.9442e+254*s^2*ws^5 + 9.7024e+262*s^2*ws^4 + 7.8861e+261*s^2*ws^3 + 2.5095e+270*s^2*ws^2 + 4.8895e+269*s^2*ws + 1.9905e+278*s^2 - 3.3148e+220*s*ws^12 + 2.6899e+223*s*ws^11 + 1.0516e+232*s*ws^10 + 2.0331e+237*s*ws^9 - 1.0373e+244*s*ws^8 + 6.8464e+248*s*ws^7 + 3.4325e+255*s*ws^6 - 6.7804e+256*s*ws^5 + 1.4597e+265*s*ws^4 - 4.3488e+264*s*ws^3 - 2.3319e+273*s*ws^2 + 7.6327e+271*s*ws + 3.1084e+280*s - 2.6732e+228*ws^12 + 2.1693e+231*ws^11 + 2.6488e+240*ws^10 - 1.4429e+243*ws^9 - 8.767e+251*ws^8 + 2.4141e+254*ws^7 + 9.6974e+262*ws^6 - 3.6803e+262*ws^5 - 1.497e+271*ws^4 + 4.8818e+269*ws^3 + 1.9883e+278*ws^2 + 2.9811e+273*ws + 1.2142e+282))/((1.9723e+81*ws^6 + 1.6586e+70*ws^5 - 1.3024e+93*ws^4 + 1.9498e+95*ws^3 + 2.1636e+104*ws^2 + 7.694e+89*ws - 3.024e+112)*(1.8535e+158*s^15 + 2.3223e+162*s^14 + 9.2675e+158*s^13*ws^2 + 1.8385e+170*s^13 + 1.0063e+163*s^12*ws^2 + 1.5501e+174*s^12 + 1.8535e+159*s^11*ws^4 + 5.5157e+170*s^11*ws^2 - 4.8956e+152*s^11*ws + 6.0797e+181*s^11 + 1.703e+163*s^10*ws^4 + 4.6214e+174*s^10*ws^2 - 8.0497e+160*s^10*ws + 2.6358e+185*s^10 + 1.8535e+159*s^9*ws^6 + 6.1289e+170*s^9*ws^4 - 9.7911e+152*s^9*ws^3 + 1.418e+182*s^9*ws^2 - 5.0127e+164*s^9*ws + 6.7076e+192*s^9 + 1.3934e+163*s^8*ws^6 + 4.607e+174*s^8*ws^4 - 4.0254e+160*s^8*ws^3 + 7.7101e+185*s^8*ws^2 - 2.0666e+183*s^8*ws + 1.6298e+195*s^8 + 9.2675e+158*s^7*ws^8 + 3.6774e+170*s^7*ws^6 + 6.0799e+181*s^7*ws^4 - 1.1547e+165*s^7*ws^3 + 2.0084e+193*s^7*ws^2 - 7.7771e+185*s^7*ws + 1.5541e+201*s^7 + 5.4188e+162*s^6*ws^8 + 1.5552e+174*s^6*ws^6 + 2.0123e+161*s^6*ws^5 + 7.4863e+185*s^6*ws^4 - 5.9779e+183*s^6*ws^3 + 3.717e+195*s^6*ws^2 - 3.6532e+193*s^6*ws + 3.677e+203*s^6 + 1.8535e+158*s^5*ws^10 + 1.8385e+170*s^5*ws^8 + 9.7911e+152*s^5*ws^7 - 6.07e+181*s^5*ws^6 - 8.0586e+164*s^5*ws^5 + 2.0074e+193*s^5*ws^4 - 1.0659e+186*s^5*ws^3 + 1.5545e+201*s^5*ws^2 - 5.7145e+195*s^5*ws + 9.4011e+208*s^5 + 7.7411e+161*s^4*ws^10 + 2.43e+172*s^4*ws^8 + 2.0124e+161*s^4*ws^7 + 2.3803e+185*s^4*ws^6 - 5.7559e+183*s^4*ws^5 + 2.6341e+195*s^4*ws^4 - 3.6377e+193*s^4*ws^3 + 1.22e+203*s^4*ws^2 - 5.3887e+200*s^4*ws + 2.2066e+211*s^4 + 6.1275e+169*s^3*ws^10 + 4.8956e+152*s^3*ws^9 - 4.0486e+181*s^3*ws^8 - 1.5239e+164*s^3*ws^7 + 6.6943e+192*s^3*ws^6 + 1.6657e+185*s^3*ws^5 + 5.2167e+200*s^3*ws^4 + 1.0076e+193*s^3*ws^3 + 1.3493e+208*s^3*ws^2 - 8.4161e+202*s^3*ws + 1.0712e+216*s^3 + 4.8461e+171*s^2*ws^10 + 4.0248e+160*s^2*ws^9 - 3.1534e+183*s^2*ws^8 - 1.8447e+183*s^2*ws^7 + 5.4181e+194*s^2*ws^6 + 3.6831e+193*s^2*ws^5 + 1.1917e+203*s^2*ws^4 - 1.0772e+201*s^2*ws^3 - 1.148e+211*s^2*ws^2 - 3.2874e+204*s^2*ws + 2.5068e+218*s^2 + 4.7489e+177*s*ws^10 + 3.795e+160*s*ws^9 - 3.1449e+189*s*ws^8 + 4.548e+185*s*ws^7 + 5.2104e+200*s*ws^6 + 5.7246e+195*s*ws^5 - 8.046e+208*s*ws^4 - 8.4161e+202*s*ws^3 + 1.0678e+216*s*ws^2 - 4.9815e+197*s*ws + 1.9579e+220*s + 3.7101e+179*ws^10 + 3.12e+168*ws^9 - 2.45e+191*ws^8 + 3.6677e+193*ws^7 + 4.0703e+202*ws^6 - 5.3837e+200*ws^5 - 6.2858e+210*ws^4 - 3.2874e+204*ws^3 + 8.3493e+217*ws^2 - 1.2972e+199*ws + 5.0985e+221))

Fig: taken from https://doi.org/10.1109/TPEL.2015.2464351

Answer 1

[Aquatris]

Here is one way:

% undamped system where w defines the natural frequency and s is the
% s-domain frequency grid variable
sys = @(s,w) 1./((s*1i).^2+2*0.05*w.*(s*1i)+w.^2); % create your tf as a function of s and w

wV = 10:1:30; % define w values you want
sV = (0:1:50);% define desired frequency 

[WW,SS] = meshgrid(wV,sV);
sysV = sys(SS,WW); % evaluate your transfer function in desired sV and wV

% plot
surf(WW,SS,abs(sysV))
xlabel('W')
ylabel('Freq')
zlabel('Magnitude')

Edit: Here is another way for your transfer function, but I think something is wrong with your derivation since the values blow up real when you evaluate:

syms ws s
G_vd =-(1.6872e+30*ws*(2.0616e+11*ws^2 - 1.5981e+19)*(- 1.0434e+209*s^14*ws^2 + 9.7644e+213*s^14*ws + 3.4484e+220*s^14 - 1.2991e+213*s^13*ws^2 + 1.2158e+218*s^13*ws + 4.2937e+224*s^13 - 5.2168e+209*s^12*ws^4 + 3.9142e+214*s^12*ws^3 + 6.8734e+220*s^12*ws^2 + 1.2884e+226*s^12*ws + 3.4203e+232*s^12 - 5.6241e+213*s^11*ws^4 + 4.8516e+218*s^11*ws^3 + 9.9349e+224*s^11*ws^2 + 9.4516e+229*s^11*ws + 2.8572e+236*s^11 - 1.0434e+210*s^10*ws^6 + 5.8925e+214*s^10*ws^5 + 3.3573e+220*s^10*ws^4 + 4.185e+226*s^10*ws^3 + 6.8393e+232*s^10*ws^2 + 3.2035e+237*s^10*ws + 1.131e+244*s^10 - 9.5051e+213*s^9*ws^6 + 7.2633e+218*s^9*ws^5 + 5.6111e+224*s^9*ws^4 + 3.0823e+230*s^9*ws^3 + 7.0634e+236*s^9*ws^2 + 9.3038e+240*s^9*ws + 4.8152e+247*s^9 - 1.0434e+210*s^8*ws^8 + 3.9566e+214*s^8*ws^7 - 1.3383e+219*s^8*ws^6 + 4.8276e+226*s^8*ws^5 + 3.4139e+232*s^8*ws^4 + 9.5864e+237*s^8*ws^3 + 2.2626e+244*s^8*ws^2 + 3.8015e+246*s^8*ws + 1.2478e+255*s^8 - 7.7621e+213*s^7*ws^8 + 4.8387e+218*s^7*ws^7 - 5.7637e+222*s^7*ws^6 + 3.5915e+230*s^7*ws^5 + 4.2122e+236*s^7*ws^4 + 2.7369e+241*s^7*ws^3 + 1.4076e+248*s^7*ws^2 + 2.6339e+249*s^7*ws + 2.0572e+257*s^7 - 5.2168e+209*s^6*ws^10 + 1.0103e+214*s^6*ws^9 - 3.5371e+220*s^6*ws^8 + 2.2567e+226*s^6*ws^7 + 3.4001e+232*s^6*ws^6 + 9.5471e+237*s^6*ws^5 + 6.907e+241*s^6*ws^4 + 9.913e+246*s^6*ws^3 + 3.7354e+255*s^6*ws^2 + 8.413e+254*s^6*ws + 2.891e+263*s^6 - 3.0095e+213*s^5*ws^10 + 1.2146e+218*s^5*ws^9 + 1.3223e+224*s^5*ws^8 + 1.7168e+230*s^5*ws^7 - 1.3395e+236*s^5*ws^6 + 2.6851e+241*s^5*ws^5 + 1.3718e+248*s^5*ws^4 + 2.9996e+249*s^5*ws^3 + 3.9942e+257*s^5*ws^2 + 2.5473e+257*s^5*ws + 4.5819e+265*s^5 - 1.0434e+209*s^4*ws^12 + 8.4668e+211*s^4*ws^11 - 6.9202e+220*s^4*ws^10 + 3.2845e+225*s^4*ws^9 + 6.818e+232*s^4*ws^8 + 3.146e+237*s^4*ws^7 - 2.2528e+244*s^4*ws^6 + 9.4169e+246*s^4*ws^5 + 3.7335e+255*s^4*ws^4 + 8.9376e+254*s^4*ws^3 + 2.891e+263*s^4*ws^2 + 4.474e+262*s^4*ws + 1.7485e+271*s^4 - 4.276e+212*s^3*ws^12 + 3.47e+215*s^3*ws^11 + 1.3494e+224*s^3*ws^10 + 2.6238e+229*s^3*ws^9 - 1.3454e+236*s^3*ws^8 + 8.7872e+240*s^3*ws^7 + 4.4564e+247*s^3*ws^6 + 1.2063e+249*s^3*ws^5 + 1.9823e+257*s^3*ws^4 + 2.6028e+256*s^3*ws^3 + 9.8906e+262*s^3*ws^2 + 8.3904e+264*s^3*ws + 2.7391e+273*s^3 - 3.4492e+220*s^2*ws^12 + 2.799e+223*s^2*ws^11 + 3.4125e+232*s^2*ws^10 - 1.8193e+235*s^2*ws^9 - 1.1279e+244*s^2*ws^8 + 3.3039e+246*s^2*ws^7 + 1.245e+255*s^2*ws^6 + 2.9442e+254*s^2*ws^5 + 9.7024e+262*s^2*ws^4 + 7.8861e+261*s^2*ws^3 + 2.5095e+270*s^2*ws^2 + 4.8895e+269*s^2*ws + 1.9905e+278*s^2 - 3.3148e+220*s*ws^12 + 2.6899e+223*s*ws^11 + 1.0516e+232*s*ws^10 + 2.0331e+237*s*ws^9 - 1.0373e+244*s*ws^8 + 6.8464e+248*s*ws^7 + 3.4325e+255*s*ws^6 - 6.7804e+256*s*ws^5 + 1.4597e+265*s*ws^4 - 4.3488e+264*s*ws^3 - 2.3319e+273*s*ws^2 + 7.6327e+271*s*ws + 3.1084e+280*s - 2.6732e+228*ws^12 + 2.1693e+231*ws^11 + 2.6488e+240*ws^10 - 1.4429e+243*ws^9 - 8.767e+251*ws^8 + 2.4141e+254*ws^7 + 9.6974e+262*ws^6 - 3.6803e+262*ws^5 - 1.497e+271*ws^4 + 4.8818e+269*ws^3 + 1.9883e+278*ws^2 + 2.9811e+273*ws + 1.2142e+282))/((1.9723e+81*ws^6 + 1.6586e+70*ws^5 - 1.3024e+93*ws^4 + 1.9498e+95*ws^3 + 2.1636e+104*ws^2 + 7.694e+89*ws - 3.024e+112)*(1.8535e+158*s^15 + 2.3223e+162*s^14 + 9.2675e+158*s^13*ws^2 + 1.8385e+170*s^13 + 1.0063e+163*s^12*ws^2 + 1.5501e+174*s^12 + 1.8535e+159*s^11*ws^4 + 5.5157e+170*s^11*ws^2 - 4.8956e+152*s^11*ws + 6.0797e+181*s^11 + 1.703e+163*s^10*ws^4 + 4.6214e+174*s^10*ws^2 - 8.0497e+160*s^10*ws + 2.6358e+185*s^10 + 1.8535e+159*s^9*ws^6 + 6.1289e+170*s^9*ws^4 - 9.7911e+152*s^9*ws^3 + 1.418e+182*s^9*ws^2 - 5.0127e+164*s^9*ws + 6.7076e+192*s^9 + 1.3934e+163*s^8*ws^6 + 4.607e+174*s^8*ws^4 - 4.0254e+160*s^8*ws^3 + 7.7101e+185*s^8*ws^2 - 2.0666e+183*s^8*ws + 1.6298e+195*s^8 + 9.2675e+158*s^7*ws^8 + 3.6774e+170*s^7*ws^6 + 6.0799e+181*s^7*ws^4 - 1.1547e+165*s^7*ws^3 + 2.0084e+193*s^7*ws^2 - 7.7771e+185*s^7*ws + 1.5541e+201*s^7 + 5.4188e+162*s^6*ws^8 + 1.5552e+174*s^6*ws^6 + 2.0123e+161*s^6*ws^5 + 7.4863e+185*s^6*ws^4 - 5.9779e+183*s^6*ws^3 + 3.717e+195*s^6*ws^2 - 3.6532e+193*s^6*ws + 3.677e+203*s^6 + 1.8535e+158*s^5*ws^10 + 1.8385e+170*s^5*ws^8 + 9.7911e+152*s^5*ws^7 - 6.07e+181*s^5*ws^6 - 8.0586e+164*s^5*ws^5 + 2.0074e+193*s^5*ws^4 - 1.0659e+186*s^5*ws^3 + 1.5545e+201*s^5*ws^2 - 5.7145e+195*s^5*ws + 9.4011e+208*s^5 + 7.7411e+161*s^4*ws^10 + 2.43e+172*s^4*ws^8 + 2.0124e+161*s^4*ws^7 + 2.3803e+185*s^4*ws^6 - 5.7559e+183*s^4*ws^5 + 2.6341e+195*s^4*ws^4 - 3.6377e+193*s^4*ws^3 + 1.22e+203*s^4*ws^2 - 5.3887e+200*s^4*ws + 2.2066e+211*s^4 + 6.1275e+169*s^3*ws^10 + 4.8956e+152*s^3*ws^9 - 4.0486e+181*s^3*ws^8 - 1.5239e+164*s^3*ws^7 + 6.6943e+192*s^3*ws^6 + 1.6657e+185*s^3*ws^5 + 5.2167e+200*s^3*ws^4 + 1.0076e+193*s^3*ws^3 + 1.3493e+208*s^3*ws^2 - 8.4161e+202*s^3*ws + 1.0712e+216*s^3 + 4.8461e+171*s^2*ws^10 + 4.0248e+160*s^2*ws^9 - 3.1534e+183*s^2*ws^8 - 1.8447e+183*s^2*ws^7 + 5.4181e+194*s^2*ws^6 + 3.6831e+193*s^2*ws^5 + 1.1917e+203*s^2*ws^4 - 1.0772e+201*s^2*ws^3 - 1.148e+211*s^2*ws^2 - 3.2874e+204*s^2*ws + 2.5068e+218*s^2 + 4.7489e+177*s*ws^10 + 3.795e+160*s*ws^9 - 3.1449e+189*s*ws^8 + 4.548e+185*s*ws^7 + 5.2104e+200*s*ws^6 + 5.7246e+195*s*ws^5 - 8.046e+208*s*ws^4 - 8.4161e+202*s*ws^3 + 1.0678e+216*s*ws^2 - 4.9815e+197*s*ws + 1.9579e+220*s + 3.7101e+179*ws^10 + 3.12e+168*ws^9 - 2.45e+191*ws^8 + 3.6677e+193*ws^7 + 4.0703e+202*ws^6 - 5.3837e+200*ws^5 - 6.2858e+210*ws^4 - 3.2874e+204*ws^3 + 8.3493e+217*ws^2 - 1.2972e+199*ws + 5.0985e+221));
G = matlabFunction(G_vd);
G(0,1) % s= 0, w =1

ans = NaN

wV = 0.6:0.05:1.6; % define w values you want
sV = logspace(0,5,1e2)*1i;% define desired frequency 
[WW,SS] = meshgrid(wV,sV);
sysV = G(SS,WW); % evaluate your transfer function in desired sV and wV
% plot
surf(WW,abs(SS),abs(sysV))
xlabel('W')
ylabel('Freq')
zlabel('Magnitude')

[Nitheesh R]

I'm not sure why you are getting real values. I'm enclosing the main body of code for your reference from where I'm deriving transfer function and plotting bode.

clc;
clear all;

s=tf('s');
fs=100e3;
ws=2*pi*fs;

V_in=80;
fi=(pi*11.2)/180;
D=fi/pi;
L_r=24.4e-6;
C_r=1.24e-7;
R_l=426.66;
R_r=0.1;
n=5.5;
C_o=30e-6;
L_i=107.5e-6;
C_i=30e-6;

A=[-(R_r/L_r) 0 0 0 ((4*sin(D*pi))/(pi*n*L_r)) ((4*cos(D*pi))/(pi*n*L_r)) -(1/L_r) 0 0 0 0 0 0 0 -(4/(pi*L_r));
   0 -(R_r/L_r) ws ((2*sin(D*pi))/(pi*n*L_r)) 0 0 0 -(1/L_r) 0 0 0 0 0 0 0;
   0 -ws -(R_r/L_r) ((2*cos(D*pi))/(pi*n*L_r)) 0 0 0 0 -(1/L_r) 0 0 0 -(2/(pi*L_r)) 0 0;
   0 -((4*sin(D*pi))/(pi*n*C_o)) -((4*cos(D*pi))/(pi*n*C_o)) -(1/(R_l*C_o)) 0 0 0 0 0 0 0 0 0 0 0;
   -((2*sin(D*pi))/(pi*n*C_o)) 0 0 0 -(1/(R_l*C_o)) ws 0 0 0 0 0 0 0 0 0;
   -((2*cos(D*pi))/(pi*n*C_o)) 0 0 0 -ws -(1/(R_l*C_o)) 0 0 0 0 0 0 0 0 0;
   (1/C_r) 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 (1/C_r) 0 0 0 0 0 0 ws 0 0 0 0 0 0;
   0 0 (1/C_r) 0 0 0 0 -ws 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 0 0 0 0 -(1/(2*L_i)) 0 -(2/(pi*L_i));
   0 0 0 0 0 0 0 0 0 0 0 ws 0 -(1/(2*L_i)) 0;
   0 0 0 0 0 0 0 0 0 0 -ws 0 (1/(pi*L_i)) 0 -(1/(2*L_i));
   -(1/(C_i)) 0 0 0 0 0 0 0 0 (1/(2*C_i)) 0 -(2/(pi*C_i)) 0 0 0;
   0 -(1/(C_i)) 0 0 0 0 0 0 0 0 (1/(2*C_i)) 0 0 0 ws;
   0 0 -(1/(C_i)) 0 0 0 0 0 0 -(1/(pi*C_i)) 0 (1/(2*C_i)) 0 -ws 0];

b=[0; 0; 0; 0; 0; 0; 0; 0; 0; (1/(L_i)); 0; 0; 0; 0; 0];

X=(-inv (A))*b*V_in;

Ir_0=X(1,1);
Ir_1R=X(2,1);
Ir_1I=X(3,1);
Vo_0=X(4,1);
Vo_1R=X(5,1);
Vo_1I=X(6,1);

f=[(((4*Vo_1R*cos(D*pi))/(n*L_r))-((4*Vo_1I*sin(D*pi))/(n*L_r)));((2*Vo_0*cos(D*pi))/(n*L_r));-((2*Vo_0*sin(D*pi))/(n*L_r));(((4*Ir_1I*sin(D*pi))/(n*C_o))-((4*Ir_1R*cos(D*pi))/(n*C_o)));-((2*Ir_0*cos(D*pi))/(n*C_o));((2*Ir_0*sin(D*pi))/(n*C_o));0;0;0;0;0;0;0;0;0];

q=[0 0 0 1 0 0 0 0 0 0 0 0 0 0 0];

I=[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 1 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 1 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 1 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 1 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 1 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 0 1 0 0 0 0 0;
   0 0 0 0 0 0 0 0 0 0 1 0 0 0 0;
   0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
   0 0 0 0 0 0 0 0 0 0 0 0 1 0 0;
   0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 1];

G_vd=q*(inv((s*I)-A))*f;
bode(G_vd)

[Aquatris]

I was getting real values with my example because I was taking the magnitude of the calculated values via the abs() function, which is plotted in the magnitude plot when you call the bode command. The phase, which can be found by the angle() command, is the 2nd plot you get with bode() function.