I have a code with several plots. How do I know from wich part of the code an specific plot comes from? (Matlab)

I have a code with several plots. How do I know from wich part of the code an specific plot comes from?

This is the example of my code:

% Living Biosystems Lab
% This code was made to optimize our hidden variables
% Optimize S1 and S2 (HF) to match S (ACDC) by comparing SV expression rates (X,Y,Z)
 
addpath("NormMI/")
 
%Setup
dt = 0.1;
Tspan = 0:dt:2500;
init_ACDC = [0.0001;0;0;0];
init_HF = [0.0001;0.0001;0.5;0.5;0.5];
options= odeset('Reltol',0.001,'Abstol',0.001,'Stats','on');
%Generate ACDC
[ts_ACDC,ys_ACDC] = ode45(@ACDC_Form_OP,Tspan,init_ACDC,options);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% Optimize signs
 
%Unless specified
W11 = 0;    
W12 = 0;
W13 = 0;
W21 = 0;
W22 = 0;
W23 = 0;
W31 = 0;   
W32 = 0;    
W33 = 0;
i = 0;
%j=0
k = 3; %degrees of freedom
for W11 = linspace(-1,1,k)
    for W22 = linspace(-1,1,k)
        for W33 = linspace(-1,1,k)
for W32 = linspace(-1,1,k) %a - effect of y on z
    for W21 = linspace(-1,1,k) %b
        for W31 = linspace(-1,1,k) %c
            for W13 = linspace(-1,1,k) %d
                
                i = i + 1;
                %j=j*2
                disp(i);
                W(:,:,i) =      [W11,    W12,    W13;
                                 W21,    W22,    W23;
                                 W31,    W32,    W33];
                           
                [ts,ys] = ode45(@(t,s)Hopf_Fold_OP(t,s,squeeze(W(:,:,i))),Tspan,init_HF);
                W(:,:,i);
                figure(3);%figure(j+1)%
                plot(ts,ys(:,3:5));
                figure(4); %figure(j+2)%
                plot3(ys(:,3),ys(:,4),ys(:,5));
                %Filter - (produces artificial innaccurarte results: That
                %is because due to the extension of the end point, it
                %creates a fake extension of the data that can influence
                %the analysis innacurately. We suggest using zeroes instead
                %- preserve the actual dynamics while allow program to
                %proceed)
                if length(ys) ~= length(Tspan)
                        ys_HF(1:length(Tspan),:,i) = [ys;zeros(length(Tspan)-length(ys),5)]; %repmat(ys(end,:),length(Tspan)-length(ys),1)];
%                         i = i - 1;
                else
                    ts_HF = ts;
                    ys_HF(:,:,i) = ys;
                end                    
            
            end
        end
    end
end
        end
    end
end
 
 
%Some weights are skippped
%% Visualize
 
figure(1);
hold on
figure(2);
hold on
for ii = 1:i
    ys = squeeze(ys_HF(:,:,ii));
    figure(1); plot(ts_HF,ys(:,3:5)) 
    figure(2); plot3(ys(:,3),ys(:,4),ys(:,5));
 
end
 
%% Visualize final sol.
 
W_mod = squeeze(W(:,:,55));
[ts_HF,ys_HF] = ode45(@(t,s)Hopf_Fold_OP(t,s,W_mod),Tspan,init_HF);
 
figure(1);
hold on
figure(2);
hold on
 
     
    figure(1);  plot(ts_HF,ys_HF(:,3:5)) 
                plot(ts_ACDC,ys_ACDC(:,2:4)) 
                    symlog(gca,'y',-1.7)
                    legend('X_HF','Y_HF','Z_HF','X_ACDC','Y_ACDC','Z_ACDC')
                    title("Comparing between Temporal Expression of the ACDC and HF Model");
                    xlabel("Time");
                    ylabel("Expression");
    
    figure(2); plot3(ys_HF(:,3),ys_HF(:,4),ys_HF(:,5));
               plot3(ys_ACDC(:,2),ys_ACDC(:,3),ys_ACDC(:,4));
                   symlog(gca,'xyz',-1.7)
                   title("Comparing between Phase Portraits of the ACDC and HF Model");
                   xlabel("X Gene");
                   ylabel("Y Gene");
                   zlabel("Z Gene");
                   
%% Measure Norm
 
%%Loss for different W
for ii = 1:i
    
   
    for sv = 1:3
        %number of the state variables (specific to the current case)
        LOSS(ii,sv) = norm((squeeze(ys_HF(:,sv+2,ii)) - ys_ACDC(:,sv+1))); %derivatives
    end
    
    LOSS_MI(ii,:) = normMI(squeeze(ys_HF(:,3:5,ii)),ys_ACDC(:,2:4),10,10);%derivatives
end
 
LOSS_log = log10(LOSS);
 
[OptLoss,OptLossIndex] = min(LOSS(:,2:4));
disp(strcat('The min loss value as per the eucledean norm is: ',num2str(OptLoss),', Thus the optimal weights are: '))
disp(W(:,:,OptLossIndex))
 
LOSS_MI_log = log10(LOSS_MI);
 
[OptLossMI,OptLossIndexMI] = min(LOSS_MI);
disp(strcat('The min loss value as per the Mutual information norm is: ',num2str(OptLossMI),', Thus the optimal weights are: '))
disp(W(:,:,OptLossIndexMI))
 
% save('Results/Weight_Optimization_Results2.mat')
 
%% Generate HeatMaps (6 for each interaction between W)
 
%pattern for W or follow
%generaate sequence then match to plot
[val,~,Trail_W13] = unique(W(1,3,:)); W13label = 'W13';
[~  ,~,Trail_W21] = unique(W(2,1,:)); W21label = 'W21';
[~  ,~,Trail_W31] = unique(W(3,1,:)); W31label = 'W31';
[~  ,~,Trail_W32] = unique(W(3,2,:)); W32label = 'W32';
Trail_W = {Trail_W13 Trail_W21 Trail_W31 Trail_W32};
Wlabel = {W13label W21label W31label W32label};
 
%Compare and HeatMap
 
i = 0;
for comp1 = 1:length(Trail_W)
    for comp2 = comp1+1:length(Trail_W)
        
        i = i + 1;
        [~,Windex(:,:,i),~] = unique([Trail_W{comp1},Trail_W{comp2}],'rows');
        Wheatmap = LOSS_log(Windex(:,:,i));
        Wheatmap = reshape(Wheatmap,k,k);
        
        figure;
        heatmap((Wheatmap)) %./min(min(Wheatmap))) - 1);
        title(strcat(Wlabel(comp1),' vs ',Wlabel(comp2)));
        
    end
end

Answer 1

[Star Strider]

The only way I can think of is to plot them with distinctive line styles or markers (or some combination of these). See the plot documentation section on: Specify Line Style, Color, and Marker.