This is robust in the event you want to find the -3dB and -6dB points on more than one peak:
fyy = linspace(0, 50, 250); % Create Data
Pyy = 0.5*exp(-(fyy-15).^2) + 0.3*exp(-(fyy-40).^2); % Create Data
[pk,loc] = findpeaks(Pyy,'Npeaks',1,'SortStr','descend');
db3c = 10^(-3/10); % Relative Magnitude At -3dB
db6c = 10^(-6/10); % Relative Magnitude At -6dB
ofst = 10;
for k1 = 1:length(pk)
varmtx = [Pyy(loc(k1)-ofst:loc(k1)); fyy(loc(k1)-ofst:loc(k1)); Pyy(loc(k1):loc(k1)+ofst); fyy(loc(k1):loc(k1)+ofst)];
dBpts(k1,1:2) = interp1(varmtx(1,:), varmtx(2,:), pk(k1)*[db6c db3c], 'linear','extrap');
dBpts(k1,3:4) = interp1(varmtx(3,:), varmtx(4,:), pk(k1)*[db6c db3c], 'linear','extrap');
end
figure(1)
plot(fyy, Pyy)
hold on
for k1 = 1:length(pk)
plot(dBpts(k1,:), pk(k1)*[db6c db3c db6c db3c], 'r+')
end
hold off
grid
for k1 = 1:length(pk)
fprintf(1, '\n\t-6dB frequencies = %.3f, %.3f\n', dBpts(k1,[1 3]))
fprintf(1, '\n\t-3dB frequencies = %.3f, %.3f\n', dBpts(k1,[2 4]))
end
