Storing different values on a matrix depending on the last value of some arrays on concatenated for loops

Hello everyone,

I am evaluating a problem that requires the resolution of a second order differential equation, which must be solved multiple times within a triple for loop. Here I leave you the code and the relevant function:

   function dx=myode(t,x,Hx,Hy,OmegaR,Omegae,gamma,alpha);
     dx=zeros(2,1);
     dx(1)=x(2);
     dx(2)=-(OmegaR^2/4)*sin(4*x(1))+2*Omegae*gamma*(Hy(t)*cos(x(1))-Hx*sin(x(1)))-...
         2*Omegae*alpha*x(2);
   end

clc
clear all
close all force
outputdir='Figures';
%% First, let's write which are the initial values of our problem
% For t=0, we have that phi=x(1)=0, and that \dot{phi}=x(2)=0
xinitial=0;
dxinitial=0;
initial_conditions=[xinitial dxinitial];
%% Involved parameters
TN=1000;
beta=0.35;
mu0=4*pi*10^(-7);
gamma=2.21*10^5;
Hx=0;
alpha0=0.001;
Omegae0=1551*gamma/mu0;
Omega4parallel0=(1645*10^(-4))*gamma/mu0;
OmegaR0=sqrt(2*Omegae0*Omega4parallel0);
%% Now, let's define the arrays involved
time=[0:0.001:90].*(10^(-12));
Temperature=[0:0.1:0.5].*TN;
pulse_duration=linspace(0,70,200).*(10^(-12));
amplitude_field=linspace(0,200,200).*((10^(-4))/mu0);
%% Now, let's propose the self-consistent method
for i=1:length(Temperature)
   magnetization=(1-(Temperature(i)/TN))^(beta); 
   alpha=(alpha0/magnetization)*(1-(Temperature(i)/(3*TN)));
   OmegaR=OmegaR0*(magnetization^(4.5)); % s^(-1)
   for j=1:length(pulse_duration)
       for k=1:length(amplitude_field)
            Hy=@(t)(t<=pulse_duration(j)).*(amplitude_field(k))+...
               (t>pulse_duration(j)).*0;
           [t,sol]=ode45(@(t,x)myode(t,x,Hx,Hy,OmegaR,Omegae0,gamma,alpha),time,initial_conditions);
           lx=cos(sol(:,1));
           ly=sin(sol(:,1));
       end
   end
end

This a priori works correctly, but obviously it takes a while to compute everything, which I have not estimated. In case it was excessive, I would already try to play with the number of elements, especially the "time" array.

What I want to do with this is the following. First of all, in each iteration, I want to be able to see which integer (-1, 0, or 1) the final element of array "lx" and array "ly" is closest to. Once this is known, I would like to know if it is possible to assign to a matrix called, for example, "matrix", the following values depending on the combination of the final values of the arrays "lx" and "ly":

• If lx(end)=1 and ly(end)=0, then matrix(j,k)=0
• If lx(end)=0 and ly(end)=1, then matrix(j,k)=pi/2
• If lx(end)=-1 and ly(end)=0, then matrix(j,k)=pi
• If lx(end)=0 and ly(end)=-1, then matrix(j,k)=3*pi/2

It is not necessary to save this matrix when each loop in "i" ends, since what I want to do for each loop in "i" is to make a phase diagram through imagesc. I mean, I would like to plot something like this:

u1=figure(1)
imagesc(pulse_duration.*(10^(12)),amplitude_field,matrix')
axis xy;
clr1=colorbar('YTickLabel',{'0','\pi/2','\pi','3\pi/2'},'YTick',[0:pi/2:3*pi/2]);
set(gca,'Ydir','normal')
t1=title('$T=%g$ K',Temperature(i),'FontSize',14,'interpreter','latex')
set(u1,'Units','Inches');
posu1=get(u1,'Position');
set(u1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[posu1(3),posu1(4)])
saveas(gcf,fullfile(outputdir,['Phase_Diagram_T=%g',Temperature(i),'_K.pdf']))

I don't know if the instruction will be correctly written to show the “Temperature” value in each iteration of “i” in both the title and the name of the saved file, let me know if it is misspelled, please. In this sense, the doubts I have regarding the graphs is (a) correctly put the title and the name of the file, (b) I do not know if by transmuting the matrix and putting ydir in the normal direction and axis xy it will be enough for the values to be conveniently labeled based on the array "amplitude_field", (c) imagine that when performing the analysis based on the last element of the arrays "lx" and "ly" no value of "matrix" is 3*pi/2, if this were the case, could its inclusion in the legend be ignored?, and (d) the ideal would be that after iteration in "Temperature", all traces of the previous phase diagram disappear (as the figure will have been saved, I don't need to keep it there), so there are no problems with u1, t1, posu1 on each iteration, right? How could that be pointed out to MatLab?

Thank you for your attention.