How to optimize one second order system for three bode plot data series?

I have three different series of measurement (bode plots: 5%, 25%, 50% of maximum amplitude of a hydraulic valve) of one nonlinear system. My objective is to find the values for gain factor(K), natural frequency(w0) and damping factor(zeta) of a linear second order system that fits best to those three bode plots. So far I managed it to fit the three series of measurement separately which gives me three different values for K,w0 and zeta. My Matlab file is attached o this request. Thank you in advance.