Well, they are supposed to be different...there’d be no point in having these multiple features if they all gave the same answer 
It’s really more about the question that each feature is trying to answer. POP gives the response of the system to a Gaussian beam. FFT and Huygens PSFs are, well, Point Spread Functions: the response of a system to a point. The FFT is faster but more approximate than the Huygens, so if they give different results when well-sampled, trust the Huygens result.
The final difference is that the FFT and Huygens PSFs only consider the diffraction is propagating from the exit pupil to the image. This is an excellent assumption, but POP will account for diffraction effects from any surface. In most imaging systems that’s not important, but it might be, and it’s also why POP has to be much higher sampled in general than the other two methods.