How to combine three different figures for three different vectors in one figure?

Here, for a = 9.1 or 3.99 or 1.1 and b = 4 or 2 or 2 i can get one single plot where i am verying another variable segma. But i want to combine them in one figure. so that i can compare them from one plot.

clc
clear variables
close all
M=16;
segma_N = 10.^-7;
p_db = 0:5:30; 
p = 10.^(p_db/10);%power
R = 0.5; %Optoelectronic conversion factor 

a = 9.1;  %a=[9.1  3.99  1.1];
b = 4;    %b=[4    2      2
];


rng(1); 
I = abs(random('normal',1,2,1,100,1))*1e-8;

%segma =[(0.0647*0.115129255)  (0.7360*0.115129255)  (4.2850*0.115129255)];
segma =[0.4472  1.264 1.870];
nz = zeros(1, length(p_db));
%ze = zeros(1, length(p_db));


for k=1:length(segma)
    
    for j = 1:length(nz) 
        ynz = qfunc(sqrt(2*((p(j).*R.*I)./segma_N).^2)).*nzpdf(I,segma(k)); %Bpsk
        nz(k,j) = trapz(I, ynz);
        
    end 
	
end

% Plot
figure
semilogy(p_db,nz(1,:),'-ks','MarkerFaceColor','k')
hold on
semilogy(p_db,nz(2,:),'-ro','MarkerFaceColor','r')
semilogy(p_db,nz(3,:),'-b^','MarkerFaceColor','b')

grid on
xlabel('P [dB]')
ylabel('ABER')
title('Intial Results')
legend('Very clear air','Haze','Light fog')
xlim([0 30])
grid on

for the function nzpdf

function G = nzpdf(I,segma)
a = 9.1;
b = 4;

c = 0.2;   %Gamma %Average optical power of classic scattering component received by off-axis eddies
%d =0.8;   %Big omega prime Average optical power of coherent contributions
p_db = 0:5:30; 
P= 10.^(p_db/10);%power

z = 2;                      %Propagation Distance
%segma = (8*0.115129255);    %Atmospheric attenuation coefficient
%segma =(4.2850*0.115129255);
s = 0.3;                    %Zero boresight 
I1 = exp(-segma*z);         %Atmospheric attenuation

omega=1.3265;
b_0=0.1079; %Average power of the coupled-to-LOS scattering component
row=0.25;  %Scattering power coupled to the LOS component
d=(omega+row.*2.*b_0+2.*sqrt(2.*b_0.*omega.*row));%Big omega prime Average optical power of coherent contributio

Aperture_radius = 0.1;      %Aperture_radius
segma_s = 0.2;              %Jitter standard deviation
omega_z =  2.5;             %Beam Width
  
v = (sqrt(pi)*Aperture_radius)/(sqrt(2)*omega_z);

omega_zeq_2 = (omega_z^2)*(sqrt(pi)*erf(v))/(2*v*exp(-(v^2))); 

g = omega_zeq_2/(2*segma_s);

A0 = (erf(v))^2;

A = (2*(a^(a/2))/((c^(1+a/2))*gamma(a)))*(((c*b)/(c*b+d))^(b+a/2));

First = (2*pi*(g^2)*A*exp((-(s^2))/(2*(segma_s^2))))/(omega_zeq_2);

count1 = 0;
for k= 1:b

    ak = (nchoosek(b-1,k-1))*(((c*b+d)^(1-k/2))/gamma(k))*((d/c)^(k-1))*((a/b)^(k/2));
    
    Summation_1 = ak*(I.^(((a+k)/2)-1)) / ( ((A0*I1)^((a+k)/2)) * sin(pi*(a-k)) );
    
    count2 = 0;
    for p = 0:length(P)
        Part1 = (((a*b.*I)/((c*b+d)*A0*I1)).^(p-(a-k)/2)) / (gamma(p-(a-k)+1)*factorial(floor(p)));
        Part2_1 = (-omega_zeq_2) / ((4*(p+k-(g^2))));
        Part2_2 = exp(((-omega_zeq_2)*(s^2)) /(8*(p+k-(g^2))*segma_s^4));
        Sum1 = Part1 .* (Part2_1 *Part2_2);
        
        Part3 = (((a*b.*I)/((c*b+d)*A0*I1)).^(p+(a-k)/2)) / (gamma(p+(a-k)+1)*factorial(floor(p)));
        Part4_1 = (-omega_zeq_2) / ((4*(p+a-(g^2))));
        Part4_2 = exp(((-omega_zeq_2)*(s^2)) / (8*(p+a-g^2)*segma_s^4));
        Min1 = Part3.*( Part4_1*Part4_2);
        
        Summation_2 = Sum1 - Min1;
        count2 = count2 + Summation_2;
    end
    
    count1 = count1 + count2.*Summation_1;

end

G = First.*count1;

end