Confusion with cut-off frequency calculation using Thevenin equivalent of an RC low pass filter
All three methods give the same answer.

The filter with the load resistor starts from -6dB, so the '-3dB point' is -9dB on that graph. Reading the first graph by eye looks consistent with just over 30Hz. Zoom into it for better accuracy.

In the second graph, for some reason you have taken the phase curve (the dotted line), which is the same height for 4 degrees shift on the right y axis and -3dB on the left y axis. You don't see the gain curve at this frequency because the y axis scale starts at -1dB (or maybe you've not plotted the gain curve). Replot your second graph with frequency going up to 100Hz, and gain axis from 0dB to -10dB. When you do, you'll see it agrees with the theoretical 31.8Hz figure you've got from calculation.