clear;
eta=0:0.01:3;
for zeta=1:-0.01:0.01
    mf=1./sqrt((1-eta.^2).^2+(2*zeta*eta).^2);
    if max(mf)>=10
        break;
    end
end
plot(eta,mf);
title(['zeta=',num2str(zeta)]);
xlabel('f/fn');ylabel('x/xst');grid('on');