function [mag phs] = sub_mybode(sys,freqVec)
  [mag,phs1] = bode(sys,freqVec*(2*pi));
  mag = squeeze(mag);  mag = 20*log10(mag);
  phs1= squeeze(phs1); abs = exp(1i*phs1*pi/180);
  phs = angle(abs)*180/pi;