I have a plant whose transfer function has all the poles in the left half of the s-plane (obtained the poles using Wolfram Alpha). However, its bode plot gives negative phase and gain margins; I obtained the bode plot from both MATLAB and Wolfram Alpha and they agree with each other. I am not able to get around this fact. How is this possible? My characteristic polynomial is
The only thing the Barkhausen criterion (which is what gain and phase margin analysis is based on) says is that in order to be oscillating a system's loop gain must be exactly 1 + 0j. In control systems we generally assume a subtraction in there someplace and turn that into the open-loop gain with a sign change must be exactly -1.
You've just discovered that the Barkhausen criterion, by itself, cannot predict stability -- it can only predict stable oscillation.
The Nyquist stability criterion is the more general test that -- if you know the number of unstable zeros in the system -- tells you whether the system is stable. I'm going to leave it to you to do the searching (a good introductory book on classical controls should have it, as does the Internet). Basically, you plot the values of the open-loop transfer function for all frequencies, and count the number of times that -1 is encircled, then compare that to the number of unstable zeros.
Personally, I prefer to start with the system in a known-stable state (found by looking at it and saying "garsh! it ain't movin'!", or by calculating the transfer function for one tuning, etc.), and then looking for gain and phase margin changes from there.