Converting MATLAB code into python code

Hello

I have written a MATLAB code, which is given below. I want to convert this code into pythod because I am not able to run this code on my laptop, which has limited Ram and very few Cores. So, i am going to ask high performance computation authority of my college to run this code. The problem is that they don't allow software which require GUI for running the codes. So, I want to convert this code into python to run it on clusters.

Please MATLAB community help me! MATLAB will always be my favourite :)

Thanks

tic 
clc; close all; clear all;     
syms p1(t) p2(t) p3(t) alpha_by_L beta_by_L gamma_by_L   epsilon_by_L 
Dp1 = diff(p1); D2p1 = diff(p1,2); Dp2 = diff(p2); D2p2 = diff(p2,2); Dp3 = diff(p3); D2p3 = diff(p3,2);
gamma_by_L_Matrix = [0.1:0.1:5];     
beta_by_L_Matrix = [0.2:0.2:10];          
alpha_by_L_Matrix = [0.1:0.1:5]; 

path = pwd ;   % mention your path 
myfolder = 'myfolder' ;   % new folder name 
folder = mkdir([path,filesep,myfolder]) ;
path  = [path,filesep,myfolder] ;

parfor (iii = 1:1:length(gamma_by_L_Matrix))
    gamma_by_L = gamma_by_L_Matrix(iii) 
figure(1)
     for jjj = 1:1:length(beta_by_L_Matrix)
         beta_by_L = beta_by_L_Matrix(jjj)
         
     for kkk = 1:1:length(alpha_by_L_Matrix)
         alpha_by_L = alpha_by_L_Matrix(kkk) 
% matrix terms 
AA = 1/2 + (alpha_by_L)*(sin(pi*beta_by_L*t))^2;
BB = (alpha_by_L)*sin(pi*beta_by_L*t)*sin(2*pi*beta_by_L*t);
CC = (alpha_by_L)*sin(pi*beta_by_L*t)*sin(3*pi*beta_by_L*t);
DD = 1/2 + (alpha_by_L)*(sin(2*pi*beta_by_L*t))^2; 
EE = (alpha_by_L)*sin(3*pi*beta_by_L*t)*sin(2*pi*beta_by_L*t);
FF = 1/2 + (alpha_by_L)*(sin(3*pi*beta_by_L*t))^2;
% matrix terms 
GG = (pi^2)/2 + (gamma_by_L)*(sin(pi*beta_by_L*t))^2;
HH = (gamma_by_L)*sin(pi*beta_by_L*t)*sin(2*pi*beta_by_L*t);  
II = (gamma_by_L)*sin(3*pi*beta_by_L*t)*sin(pi*beta_by_L*t);
JJ = (4*pi^2)/2 + (gamma_by_L)*(sin(2*pi*beta_by_L*t))^2;
KK = (gamma_by_L)*sin(2*pi*beta_by_L*t)*sin(3*pi*beta_by_L*t);
LL = (9*pi^2)/2 + (gamma_by_L)*(sin(3*pi*beta_by_L*t))^2;
% RHS    
MM = (epsilon_by_L)*sin(pi*beta_by_L*t);
NN = (epsilon_by_L)*sin(2*pi*beta_by_L*t); 
OO = (epsilon_by_L)*sin(3*pi*beta_by_L*t);
    
 Eq1 = AA*diff(p1,t,2) + BB*diff(p2,t,2) + CC*diff(p3,t,2) +  GG*p1 + HH*p2 + II*p3 == 0; % Equation 1
 Eq2 = BB*diff(p1,t,2) + DD*diff(p2,t,2) + EE*diff(p3,t,2) +  HH*p1 + JJ*p2 + KK*p3 == 0; % Equation 2
 Eq3 = CC*diff(p1,t,2) + EE*diff(p2,t,2) + FF*diff(p3,t,2) +  II*p1 + KK*p2 + LL*p3 == 0; % Equation 3

[V,S] = odeToVectorField(Eq1, Eq2, Eq3); 
ftotal = matlabFunction(V, 'Vars',{'t','Y'});           % Using readymade MATLAB function to solve using ODE 45
interval = [0 2/beta_by_L];                                  % Time Interval to run the program 

for i = 1:6 
   if i == 1 
       y0 = [1 ;0 ;0 ;0 ;0 ;0 ] ; % First IC
    
   elseif i == 2
       y0 = [0 ;1 ;0 ;0 ;0 ;0 ] ; % Second IC
 
   elseif i == 3
       y0 = [0 ;0 ;1 ;0 ;0 ;0 ] ;
   
   elseif i == 4
       y0 = [0 ;0 ;0 ;1 ;0 ;0 ] ;
   
   elseif i == 5
       y0 = [0 ;0 ;0 ;0 ;1 ;0 ] ; 
   
   elseif i == 6
       y0 = [0 ;0 ;0 ;0 ;0 ;1 ] ;     
   end
xSol = ode45( @(t,Y) ftotal(t,Y),interval,y0);        % Using ODE 45 to solve state space form of ODE
tValues = linspace(interval(1),interval(2),1000);     % dividing time interval      
x1Values = deval(xSol,tValues,1); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
x2Values = deval(xSol,tValues,2); % number 2 denotes second solution likewise you can mention 2 ,3 & 4 for the next three solutions
x3Values = deval(xSol,tValues,3); 
x4Values = deval(xSol,tValues,4);
x5Values = deval(xSol,tValues,5);
x6Values = deval(xSol,tValues,6);
if i == 1
col_11 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_12 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
 
elseif i == 2
col_21 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_22 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
  
elseif i == 3
col_31 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_32 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

elseif i == 4
col_41 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_42 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

elseif i == 5
col_51 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_52 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

elseif i == 6
col_61 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
col_62 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
 
end
end     
A1 = [col_11 col_21 col_31 col_41 col_51 col_61];
A2 = [col_12 col_22 col_32 col_42 col_52 col_62] ;
A = (inv(A1))*A2  ;
eigval = eig(A);
abseigval = abs(eigval);
iteval = gamma_by_L

if (abseigval(1)<=1.00 && abseigval(2)<=1.00 && abseigval(3)<=1.00 && abseigval(4)<=1.00 && abseigval(5)<=1.00 && abseigval(6)<=1.00)
             plot(beta_by_L,alpha_by_L,'k.', 'MarkerSize', 10)
             xlabel('beta_by_L')
             ylabel('alpha_by_L')
title("gamma_by_L " + gamma_by_L + " ")

            hold on
else       
            plot(beta_by_L,alpha_by_L,'r.', 'MarkerSize', 10)
            xlabel('beta_by_L')
             ylabel('alpha_by_L')
title("gamma_by_L " + gamma_by_L + " ")

            hold on          
    
end 
     end 
     end
     
     hold off
    xlim([0 10])
    ylim([0 5])
    iii         

    
end
    toc 
       

[Mike Croucher]

You don't need a GUI to run MATLAB code. You can run it in batch mode just like you can Python. A huge part of my job is to work with researchers who want to run MATLAB effectively on High Performance Computing clusters.

Your first issue with this code is 'What results do you want to save?' I'm guessing the variables beta_by_L,alpha_by_L since you are plotting them?

[aakash dewangan]

Thank you very much Mike for reply :)

I am glad to know that I can use MATLAB in high performance computation clusters :)

From this code, I will be getting plots/graphs of "beta_by_L vs alpha_by_L" for different values of "gamma_by_L". For different values of gamma_by_L, I will get different graphs. I want all those graphs in my result.

Please help me and guide me to run this on high performance computing cluster.

I really need your help.

Answer 1

[Mike Croucher]

Hi Aakash

To get this running on a HPC system, it would be beneficial to know the details of the HPC system. As I said in my comment to your question, you can run MATLAB in batch mode (that is without a GUI). The relevant section of the documentation is Start MATLAB program from Linux system prompt - MATLAB - MathWorks United Kingdom and the command switch you want is matlab -batch yourcode.m

Before going further down that route, however, I suggest some changes to your code.

Note that I have only done limited testing here and did the smallest amount of changes I could to get the speed-ups. You should check that it gives the results you need!

There's more that could be done to improve the readability of this code that might not affect the speed. I deliberately didn't do any of that.

On my machine, the first iteration of iii went down from 2942 to 20 seconds. A speed-up of 147 times.

The full run took 227 seconds to produce 50 plots on an 8 core system.

The only 'test' I've done is to check that the first plot given corresponding to iii=1 is the same in both your original code and the new one....and this was done by eye. Please check yourself that you get the results you expect.

The original code gave the following in 2942 seconds for iii=1

the final version of the code gave the following in 20 seconds for iii=1

The first idea is to move all of the slow symbolic calculations outside of the main loop. This provides most of the speed-up, around 70x or so

%% new version of code that speeds up the symbolic section
%% this is not the fastest version so far
tic
clc; close all; clear all;
syms p1(t) p2(t) p3(t) alpha_by_L beta_by_L gamma_by_L   epsilon_by_L
syms Sgamma_by_L Sbeta_by_L Salpha_by_L
Dp1 = diff(p1); D2p1 = diff(p1,2); Dp2 = diff(p2); D2p2 = diff(p2,2); Dp3 = diff(p3); D2p3 = diff(p3,2);
gamma_by_L_Matrix = [0.1:0.1:5];
beta_by_L_Matrix = [0.2:0.2:10];
alpha_by_L_Matrix = [0.1:0.1:5];

%%You created an output folder here but you are not saving
%%anything to it.  I've chaged it a little and save the figure of each
%%iteration.
resultsfolder = "myfolder";
mkdir(resultsfolder);

% Create the symbolic matrix, and a slightly modified ftotal once instead
% of every iteration.  The new ftotal function can be used every iteration
% without being created from scratch
diff_p1_t_2 = diff(p1,t,2);
diff_p2_t_2 = diff(p2,t,2);
diff_p3_t_2 = diff(p3,t,2);
AA = 1/2 + (Salpha_by_L)*(sin(pi*Sbeta_by_L*t))^2;
BB = (Salpha_by_L)*sin(pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
CC = (Salpha_by_L)*sin(pi*Sbeta_by_L*t)*sin(3*pi*Sbeta_by_L*t);
DD = 1/2 + (Salpha_by_L)*(sin(2*pi*Sbeta_by_L*t))^2;
EE = (Salpha_by_L)*sin(3*pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
FF = 1/2 + (Salpha_by_L)*(sin(3*pi*Sbeta_by_L*t))^2;
% matrix terms
GG = (pi^2)/2 + (Sgamma_by_L)*(sin(pi*Sbeta_by_L*t))^2;
HH = (Sgamma_by_L)*sin(pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
II = (Sgamma_by_L)*sin(3*pi*Sbeta_by_L*t)*sin(pi*Sbeta_by_L*t);
JJ = (4*pi^2)/2 + (Sgamma_by_L)*(sin(2*pi*Sbeta_by_L*t))^2;
KK = (Sgamma_by_L)*sin(2*pi*Sbeta_by_L*t)*sin(3*pi*Sbeta_by_L*t);
LL = (9*pi^2)/2 + (Sgamma_by_L)*(sin(3*pi*Sbeta_by_L*t))^2;
% RHS
MM = (epsilon_by_L)*sin(pi*Sbeta_by_L*t);
NN = (epsilon_by_L)*sin(2*pi*Sbeta_by_L*t);
OO = (epsilon_by_L)*sin(3*pi*Sbeta_by_L*t);

Eq1 = AA*diff_p1_t_2 + BB*diff_p2_t_2 + CC*diff_p3_t_2 +  GG*p1 + HH*p2 + II*p3 == 0; % Equation 1
Eq2 = BB*diff_p1_t_2 + DD*diff_p2_t_2 + EE*diff_p3_t_2 +  HH*p1 + JJ*p2 + KK*p3 == 0; % Equation 2
Eq3 = CC*diff_p1_t_2 + EE*diff_p2_t_2 + FF*diff_p3_t_2 +  II*p1 + KK*p2 + LL*p3 == 0; % Equation 3
[V,S] = odeToVectorField(Eq1, Eq2, Eq3);
ftotal = matlabFunction(V, 'Vars',{'t','Y','Salpha_by_L','Sbeta_by_L','Sgamma_by_L'});           % Using readymade MATLAB function to solve using ODE 45

parfor iii = 1:length(gamma_by_L_Matrix)
    gamma_by_L = gamma_by_L_Matrix(iii);
    myfig = figure(1)
    for jjj = 1:1:length(beta_by_L_Matrix)
        beta_by_L = beta_by_L_Matrix(jjj);
        for kkk = 1:1:length(alpha_by_L_Matrix)
            alpha_by_L = alpha_by_L_Matrix(kkk);
            interval = [0 2/beta_by_L];                                  % Time Interval to run the program
            %% This section could be rewritten a lot but I leave it the same
            %  to make it easier for you to follow the changes done so far
            for i = 1:6
                if i == 1
                    y0 = [1 ;0 ;0 ;0 ;0 ;0 ]; % First IC
                elseif i == 2
                    y0 = [0 ;1 ;0 ;0 ;0 ;0 ]; % Second IC
                elseif i == 3
                    y0 = [0 ;0 ;1 ;0 ;0 ;0 ];
                elseif i == 4
                    y0 = [0 ;0 ;0 ;1 ;0 ;0 ];
                elseif i == 5
                    y0 = [0 ;0 ;0 ;0 ;1 ;0 ];
                elseif i == 6
                    y0 = [0 ;0 ;0 ;0 ;0 ;1 ];
                end
                xSol = ode45( @(t,Y) ftotal(t,Y,alpha_by_L,beta_by_L,gamma_by_L),interval,y0);        % Using ODE 45 to solve state space form of ODE
                tValues = linspace(interval(1),interval(2),1000);     % dividing time interval
                x1Values = deval(xSol,tValues,1); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
                x2Values = deval(xSol,tValues,2); % number 2 denotes second solution likewise you can mention 2 ,3 & 4 for the next three solutions
                x3Values = deval(xSol,tValues,3);
                x4Values = deval(xSol,tValues,4);
                x5Values = deval(xSol,tValues,5);
                x6Values = deval(xSol,tValues,6);
                if i == 1
                    col_11 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
                    col_12 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
                elseif i == 2
                    col_21 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
                    col_22 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

                elseif i == 3
                    col_31 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
                    col_32 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
                elseif i == 4
                    col_41 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
                    col_42 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
                elseif i == 5
                    col_51 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
                    col_52 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
                elseif i == 6
                    col_61 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
                    col_62 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
                end
            end
            A1 = [col_11 col_21 col_31 col_41 col_51 col_61];
            A2 = [col_12 col_22 col_32 col_42 col_52 col_62] ;
            A = (inv(A1))*A2  ;
            eigval = eig(A);
            abseigval = abs(eigval);
            iteval = gamma_by_L;
            %% You make a plot here for every iteration of the iii loop but are not saving it.  
            %% Is it these plots that you want to save?
            if (abseigval(1)<=1.00 && abseigval(2)<=1.00 && abseigval(3)<=1.00 && abseigval(4)<=1.00 && abseigval(5)<=1.00 && abseigval(6)<=1.00)
                plot(beta_by_L,alpha_by_L,'k.', 'MarkerSize', 10)
                xlabel('beta_by_L')
                ylabel('alpha_by_L')
                title("gamma_by_L " + gamma_by_L + " ")
                hold on
            else
                plot(beta_by_L,alpha_by_L,'r.', 'MarkerSize', 10)
                xlabel('beta_by_L')
                ylabel('alpha_by_L')
                title("gamma_by_L " + gamma_by_L + " ")
                hold on
            end
        end
    end
    hold off
    %% An example of how to save the plots to a unique filename per iteration
    filename = resultsfolder+filesep+"result_"+iii+".fig";
    savefig(filename);
    %% I would suggest saving a numerical result as well per iteration.  Have a think about what that should be
    xlim([0 10])
    ylim([0 5])
    iii
end
toc

The next thing you can do is vectorise your calls to deval. i.e. change

% original deval section
x1Values = deval(xSol,tValues,1); % number 1 denotes first solution likewise you can mention 2 ,3 & 4 for the next three solutions
x2Values = deval(xSol,tValues,2); % number 2 denotes second solution likewise you can mention 2 ,3 & 4 for the next three solutions
x3Values = deval(xSol,tValues,3);
x4Values = deval(xSol,tValues,4);
x5Values = deval(xSol,tValues,5);
x6Values = deval(xSol,tValues,6);
if i == 1
    col_11 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
    col_12 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

elseif i == 2
    col_21 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
    col_22 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];

elseif i == 3
    col_31 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
    col_32 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 4
    col_41 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
    col_42 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 5
    col_51 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
    col_52 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
elseif i == 6
    col_61 = [x1Values(1);x2Values(1);x3Values(1);x4Values(1);x5Values(1);x6Values(1)];
    col_62 = [x1Values(end);x2Values(end);x3Values(end);x4Values(end);x5Values(end);x6Values(end)];
end

to

%% New deval section
xvalues = deval(xSol,tValues,1:6);
 
if i == 1
    col_11 = xvalues(:,1);
    col_12 = xvalues(:,end);
elseif i == 2
    col_21 = xvalues(:,1);
    col_22 = xvalues(:,end);
elseif i == 3
    col_31 = xvalues(:,1);
    col_32 = xvalues(:,end);
elseif i == 4
    col_41 = xvalues(:,1);
    col_42 = xvalues(:,end);
elseif i == 5
    col_51 = xvalues(:,1);
    col_52 =xvalues(:,end);
elseif i == 6
    col_61 =  xvalues(:,1);
    col_62 = xvalues(:,end);
end

This gave another 2x speed-up on my machine. One iteration of iii took 20 seconds now with the full calculation of 50 iterations taking 227 seconds on my 8 core machine. The full version of the new code is:

%% The fastest version so far
tic
clc; close all; clear all;
syms p1(t) p2(t) p3(t) alpha_by_L beta_by_L gamma_by_L   epsilon_by_L
syms Sgamma_by_L Sbeta_by_L Salpha_by_L
Dp1 = diff(p1); D2p1 = diff(p1,2); Dp2 = diff(p2); D2p2 = diff(p2,2); Dp3 = diff(p3); D2p3 = diff(p3,2);
gamma_by_L_Matrix = [0.1:0.1:5];
beta_by_L_Matrix = [0.2:0.2:10];
alpha_by_L_Matrix = [0.1:0.1:5];

%%You created an output folder here but you are not saving
%%anything to it.  I've chaged it a little and save the figure of each
%%iteration.
resultsfolder = "myfolder";
mkdir(resultsfolder);

% Create the symbolic matrix, and a slightly modified ftotal once instead
% of every iteration.  The new ftotal function can be used every iteration
% without being created from scratch
diff_p1_t_2 = diff(p1,t,2);
diff_p2_t_2 = diff(p2,t,2);
diff_p3_t_2 = diff(p3,t,2);
AA = 1/2 + (Salpha_by_L)*(sin(pi*Sbeta_by_L*t))^2;
BB = (Salpha_by_L)*sin(pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
CC = (Salpha_by_L)*sin(pi*Sbeta_by_L*t)*sin(3*pi*Sbeta_by_L*t);
DD = 1/2 + (Salpha_by_L)*(sin(2*pi*Sbeta_by_L*t))^2;
EE = (Salpha_by_L)*sin(3*pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
FF = 1/2 + (Salpha_by_L)*(sin(3*pi*Sbeta_by_L*t))^2;
% matrix terms
GG = (pi^2)/2 + (Sgamma_by_L)*(sin(pi*Sbeta_by_L*t))^2;
HH = (Sgamma_by_L)*sin(pi*Sbeta_by_L*t)*sin(2*pi*Sbeta_by_L*t);
II = (Sgamma_by_L)*sin(3*pi*Sbeta_by_L*t)*sin(pi*Sbeta_by_L*t);
JJ = (4*pi^2)/2 + (Sgamma_by_L)*(sin(2*pi*Sbeta_by_L*t))^2;
KK = (Sgamma_by_L)*sin(2*pi*Sbeta_by_L*t)*sin(3*pi*Sbeta_by_L*t);
LL = (9*pi^2)/2 + (Sgamma_by_L)*(sin(3*pi*Sbeta_by_L*t))^2;
% RHS
MM = (epsilon_by_L)*sin(pi*Sbeta_by_L*t);
NN = (epsilon_by_L)*sin(2*pi*Sbeta_by_L*t);
OO = (epsilon_by_L)*sin(3*pi*Sbeta_by_L*t);

Eq1 = AA*diff_p1_t_2 + BB*diff_p2_t_2 + CC*diff_p3_t_2 +  GG*p1 + HH*p2 + II*p3 == 0; % Equation 1
Eq2 = BB*diff_p1_t_2 + DD*diff_p2_t_2 + EE*diff_p3_t_2 +  HH*p1 + JJ*p2 + KK*p3 == 0; % Equation 2
Eq3 = CC*diff_p1_t_2 + EE*diff_p2_t_2 + FF*diff_p3_t_2 +  II*p1 + KK*p2 + LL*p3 == 0; % Equation 3
[V,S] = odeToVectorField(Eq1, Eq2, Eq3);
ftotal = matlabFunction(V, 'Vars',{'t','Y','Salpha_by_L','Sbeta_by_L','Sgamma_by_L'});           % Using readymade MATLAB function to solve using ODE 45

parfor iii = 1:length(gamma_by_L_Matrix)
    gamma_by_L = gamma_by_L_Matrix(iii);
    myfig = figure(1)
    for jjj = 1:1:length(beta_by_L_Matrix)
        beta_by_L = beta_by_L_Matrix(jjj);
        for kkk = 1:1:length(alpha_by_L_Matrix)
            alpha_by_L = alpha_by_L_Matrix(kkk);
            interval = [0 2/beta_by_L];                                  % Time Interval to run the program
            %% This section could be rewritten a lot but I leave it the same to make it easier to follow the changes done so far
            for i = 1:6
                if i == 1
                    y0 = [1 ;0 ;0 ;0 ;0 ;0 ]; % First IC
                elseif i == 2
                    y0 = [0 ;1 ;0 ;0 ;0 ;0 ]; % Second IC
                elseif i == 3
                    y0 = [0 ;0 ;1 ;0 ;0 ;0 ];
                elseif i == 4
                    y0 = [0 ;0 ;0 ;1 ;0 ;0 ];
                elseif i == 5
                    y0 = [0 ;0 ;0 ;0 ;1 ;0 ];
                elseif i == 6
                    y0 = [0 ;0 ;0 ;0 ;0 ;1 ];
                end
                xSol = ode45( @(t,Y) ftotal(t,Y,alpha_by_L,beta_by_L,gamma_by_L),interval,y0);        % Using ODE 45 to solve state space form of ODE
                tValues = linspace(interval(1),interval(2),1000);     % dividing time interval
                xvalues = deval(xSol,tValues,1:6);
                 
                if i == 1
                    col_11 = xvalues(:,1);
                    col_12 = xvalues(:,end);
                elseif i == 2
                    col_21 = xvalues(:,1);
                    col_22 = xvalues(:,end);
                elseif i == 3
                    col_31 = xvalues(:,1);
                    col_32 = xvalues(:,end);
                elseif i == 4
                    col_41 = xvalues(:,1);
                    col_42 = xvalues(:,end);
                elseif i == 5
                    col_51 = xvalues(:,1);
                    col_52 =xvalues(:,end);
                elseif i == 6
                    col_61 =  xvalues(:,1);
                    col_62 = xvalues(:,end);
                end
            end
            A1 = [col_11 col_21 col_31 col_41 col_51 col_61];
            A2 = [col_12 col_22 col_32 col_42 col_52 col_62] ;
            A = (inv(A1))*A2  ;
            eigval = eig(A);
            abseigval = abs(eigval);
            iteval = gamma_by_L;
            %% You make a plot here for every iteration of the iii loop but are not saving it.  
            %% Is it these plots that you want to save?
            if (abseigval(1)<=1.00 && abseigval(2)<=1.00 && abseigval(3)<=1.00 && abseigval(4)<=1.00 && abseigval(5)<=1.00 && abseigval(6)<=1.00)
                plot(beta_by_L,alpha_by_L,'k.', 'MarkerSize', 10)
                xlabel('beta_by_L')
                ylabel('alpha_by_L')
                title("gamma_by_L " + gamma_by_L + " ")
                hold on
            else
                plot(beta_by_L,alpha_by_L,'r.', 'MarkerSize', 10)
                xlabel('beta_by_L')
                ylabel('alpha_by_L')
                title("gamma_by_L " + gamma_by_L + " ")
                hold on
            end
        end
    end
    hold off
    %% An example of how to save the plots to a unique filename per iteration
    filename = resultsfolder+filesep+"result_"+iii+".fig";
    savefig(filename);
    %% I would suggest saving a numerical result as well per iteration.  Have a think about what that should be
    xlim([0 10])
    ylim([0 5])
    iii
end
toc

Hope this helps,

Mike

[aakash dewangan]

O My GOD!!!! O My GOD!!!! O My GOD!!!! O My GOD!!!! O My GOD!!!! O My GOD!!!! O My GOD!!!!

I can't express my feeling in words. You are legend in the MATLAB sir. Thank you so much for your help. thank you very very very very very much sir. It works exactly in the same way, and the results are also same. I don't know how to thank you for this help. You have no idea what you have done for me. MATLAB is best...... And you are the best in solving the problems of MATLAB.

Best Regards,