I think you might find this answer helpful: https://www.mathworks.com/matlabcentral/answers/477222-copy-a-figure-to-a-subplot-including-all-elements#answer_388666
I think the issue is that, in your code, ax always refers to the same axes, not a copy of the template axes you created. Because of that, you keep moving the axes around, with it ending in the final spot with all plots in it.
Instead, you want to create a copy of the axes object in each subplot location, which is what the linked answer does. It's not quite the elegant solution you might expect or want, but it works.
load('data.mat');
figure(1)
ax = axes();
hold(ax);
ax.YAxis.Scale = 'log';
ax.XAxis.Scale = 'log';
%xlabel(ax, 'Frequency (Hz)', 'Interpreter', 'latex', 'FontSize', 14,'FontWeight','Bold');
%ylabel(ax, 'Power spectral density', 'Interpreter', 'latex', 'FontSize', 14,'FontWeight','Bold');
grid on
% setup top axis
ax_top = axes(); % axis to appear at top
hold(ax_top);
ax_top.XAxis.Scale = 'log';
ax_top.XAxisLocation = 'top';
ax_top.YAxisLocation = "right";
ax_top.YTick = [];
% ax_top.XDir = 'reverse';
ax_top.Color = 'none';
%xlabel(ax_top, 'Strouhal number', 'Interpreter', 'latex', 'FontSize', 14,'FontWeight','Bold');
% linking axis
linkprop([ax, ax_top],{'Units','Position','ActivePositionProperty'});
%ax.Position(1,1) = ax.Position(1);
v = 4.687276290832519e+02*12; % inches/sec
c = 16; % inches
f=[0.25 0.5 1 2 4 6 10 16 30 50 100 200 500];
st = f*c./v; % Strouhal number
% configure limits of bottom axis
ax.XLim = [f(1) f(end)];
ax.XTick = f;
ax.XAxis.MinorTick = 'on';
% configure limits and labels of top axis
ax_top.XLim = [st(1) st(end)];
% Updated code
fig2 = figure(2);
ax1 = subplot(2,2,1,'parent',fig2);
ax2 = subplot(2,2,2,'parent',fig2);
ax3 = subplot(2,2,3,'parent',fig2);
ax4 = subplot(2,2,4,'parent',fig2);
axcp1 = copyobj(ax, fig2);
set(axcp1,'Position',get(ax1,'position'));
delete(ax1);
plot(axcp1,fCompnnts(1:end/2),psd3(1:end/2))
title(axcp1,'$x/c=0.448$','Interpreter','latex')
axcp2 = copyobj(ax, fig2);
set(axcp2,'Position',get(ax2,'position'));
delete(ax2);
plot(axcp2,fCompnnts(1:(end)/2),psd4(1:end/2))
title(axcp2,'$x/c=0.498$','Interpreter','latex')
axcp3 = copyobj(ax, fig2);
set(axcp3,'Position',get(ax3,'position'));
delete(ax3);
plot(axcp3,fCompnnts(1:end/2),psd5(1:end/2))
title(axcp3,'$x/c=0.542$','Interpreter','latex')
axcp4 = copyobj(ax, fig2);
set(axcp4,'Position',get(ax4,'position'));
delete(ax4);
plot(axcp4,fCompnnts(1:end/2),psd8(1:end/2),'-h','Color',[0.6350 0.0780 0.1840])
title(axcp4,'$x/c=0.950$','Interpreter','latex')
I get an error when trying to run that code here, but it works in R2022a at least. Here is the result
