plot 95% confidence interval for a quantile-quantile plot

I'm trying to plot 95% confidence interval for the plot of two sets of simulated quantiles as part of a q-q analysis. I know how to compute the 95% confidence bounds for each set of quantiles (i.e. x quantiles or y quantiles), but what is the best way to compute the 95% confidence interval of both x and y when plotted together in a q-q plot? I don't think plotting the upper (lower) bounds of x vs. the lower (upper) bounds of y is the correct way to go about this (see my code below). I am not looking for a ci based on a line of best fit. Instead I'm looking for some way to compute the upper and lower bounds of y at each value of x or similar.

ci = 0.1;

cih1 = ceil((mc - mc.*(ci/2))); 
cih2 = floor((mc - mc.*(ci/2)));
cil1 = ceil((mc.*(ci/2)));
cil2 = floor((mc.*(ci/2)));

parfor t=1:q
    slgm(:,t) = sort(qlgm_sim(:,t));
end

qlgm_h1 = slgm(cih1,:);
qlgm_h2 = slgm(cih2,:);
qlgm_l1 = slgm(cil1,:);
qlgm_l2 = slgm(cil2,:);

qlgmh = vertcat(qlgm_h1,qlgm_h2);
hilgm= mean(qlgmh);
qlgml = vertcat(qlgm_l1,qlgm_l2);
lolgm = mean(qlgml);