Solve a second-order differential equation with constant parameters changing through an external parameter

Hello everyone,

I know that this a very vague question, but I am quite lost in this aspect, so I am asking for some help to begin to try to solve a second order differential equation. My second-order differential equation has the following functional form: d^2 phi (t)/dx^2+OmegaR^2/4*sin(4*phi(t))-2*Omegae*gamma*(Hy*cos(phi(t)-Hx*sin(phi(t)))+2*Omegae*alpha*dphi(t)/dt. I know which are my initial conditions, phi(t=0)=0, dphi(t)/dt (t=0)=0. Even I know which is the value of my coefficients OmegaR, Omegae, gamma, Hy, Hx and alpha for each value of an external temperature, let's say, temperature, Temperature variable, not involved in the second-order differential equation. So at the end, I want to solve that second-order differential equation in a certain time (t) interval, and plot it. It would be great if in a future I can do implement all the solutions for all the temperatures in a single plot indicating in a legend for which temperature each line corresponds, but that it is other stage of the problem. I have read the documentation of MatLab, and some posts, but I am quite lost, certainly. Someone can give an example and a strategy from which I could begin?

Also I would want to give a time-dependence to H_x and H_y, like being linear in a certain time interval H_x=k*t, given a k value, and after reaching a certain value, to be constant the rest of the time imposed in the process of solving the second-order differential equation.

Thank for your attention.

Edit: Here I present my approach until now

First, I create a function file, with all the parameters defined for a certain value of the external parameter temperature (Temperature), equal to 0. First, I am not sure if I have to upload the values in my function file (I mean, OmegaR, Omegae, gamma, Hy, Hx and alpha, because I the generalized case I will have to load a file which will contain these values for different value of temperature, which I have to remember that it is distinct to phi and time (t)):

   function dx=fx(t,x);
     % Global variables: Hx, Hy, OmegaR, Omegae, gamma, alpha
     % Here, x(1)=phi, x(2)=dot{phi}
     dx=zeros(2,1);
     dx(1)=x(2);
     mu0=4*pi*10^(-7); %H/m
     gamma=2.21*10^5; %m/(A*s)
     kB=1.38064852*10^(-23); %J/K
     a0=3.328*10^(-10); %m
     c=8.539*10^(-10); %m
     V=c*(a0^2); %m^3
     Hx=0;
     Hy=100*79.5774715459; %A/m
     alpha=0.001;
     muB=9.27400994*10^(-24); %(T^2*m^3)/J
     mu=4*muB; %(T^2*m^3)/J
     Ms=mu/V; %T^2/J
     Omegae=(2*gamma*abs(4*(-396*kB/V)-532*kB/V))/(mu0*Ms); %s^-1
     K4parallel=1.8548*(10^(-25))/V; %J/m^3
     Omega4parallel=(2*gamma*K4parallel)/(mu0*Ms); %s^-1
     OmegaR=sqrt(2*Omegae*Omega4parallel);
     dx(2)=-(OmegaR.^2/4).*sin(4*x(1))-2*Omegae*gamma*(Hy*cos(x(1))-Hx*sin(x(1)))+...
         2*Omegae*alpha*x(2);
   end

After that, I can solve my equation for a certain time interval:

clc
clear all
close all force
output='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];
%% Now, let's create the array of times on which we are interested
time=[0:0.001:40].*(10^(-12)); %s
%% Let's solve the second-order differential equation, where the X variable will have to columns
% Probably, the first column corresponds to phi=x(1) and the second one to
% \dot{phi}=x(2)
[T,X]=ode45(@fx,time,initial_conditions);

No error occurs. But I am not sure why X is an array of 40001x2 elements, what it is in each column? The first one phi and the second one its first derivative? Moreover, as I have said, this for the more simplified case of a single set of parameters for a certain value of the external parameter temperature, so I was wonder how to confront the situation in which I have multiple values for the aforementioned parameters.

Moreover, I am interested in how to implement the solution of my second-order differential equation without specifying a priori the values of Hx and Hy, because I want to take control on each of them each time that I run the program. It would be great if I run the program and MatLab ask me which value of Hx and Hy I want to use.

[darova]

What is the problem? ode45 and go?

[Richard Wood]

Dear Darova,

thank you for your answer. Please, see above my editted post.

Answer 1

[Star Strider]

It appears that what you want to do is to pass your data to your ODE function ‘fx’. See the documentation section on Passing Extra Parameters to understand how to do that.

Please do not use global variables!

function dx=fx(t,x);
     % Global variables: Hx, Hy, OmegaR, Omegae, gamma, alpha

Do this instead:

     function dx = fx(t, x, Hx, Hy, OmegaR, Omegae, gamma, alpha);

and your ode45 call then becomes:

         [T,X] = ode45(@(t,x)fx(t, x, Hx, Hy, OmegaR, Omegae, gamma, alpha), time, initial_conditions);

I assume that the additional variables are scalars that do not change during the integration, and that you want to run several integrations, each with different additional parameters. You can do this with a for loop, saving ‘T’ and ‘X’ as cell arrays in each iteration (easiest), then plotting them individually at the end (likely also with a loop).

[Richard Wood]

Hey! Thank you very much! That was a nice answer. Nevertheless, I have some other questions related to that that you have told me. Your assumption is correct, the variables Hx, Hy, OmegaR, Omegae, gamma and alpha are scalar for each temperature value (so they are a constant in the loop to solve the second-order differential equation). With that in mind, I can find the value of certain new variables in each temperature iteration, that it is, lx, ly, and mz. After that, I want to plot, for example, the lx variable, once I have the value of lx for each time step for each temperature value. Moreover, I want to place in a legend which color corresponds for which temperature. After that it is complete, I want to save the final plot. Obviously the attempt that I write to you it is not correct, but some initial scheme. Give me your opinion if possible:

   function dx=fx(t,x,Hx,Hy,OmegaR,Omegae,gamma,alpha);
     % Here, x represents the in-plane angle, phi, and t represents the
     % time variable
     % Global variables: Hx, Hy, OmegaR, Omegae, gamma, alpha
     % Here, x(1)=phi, x(2)=\dot{phi}
     dx=zeros(2,1);
     dx(1)=x(2);
     dx(2)=-(OmegaR.^2/4).*sin(4*x(1))-2*Omegae*gamma*(Hy*cos(x(1))-Hx*sin(x(1)))+...
         2*Omegae*alpha*x(2);
   end

   clc
clear all
close all force
output='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];
% Also, I need to define a parameter, Omegae
mu0=4*pi*10^(-7); %H/m
gamma=2.21*10^5; %m/(A*s)
kB=1.38064852*10^(-23); %J/K
a0=3.328*10^(-10); %m
c=8.539*10^(-10); %m
V=c*(a0^2); %m^3
Hx=0;
Hy=100*79.5774715459; %A/m
alpha=0.001;
muB=9.27400994*10^(-24); %(T^2*m^3)/J
mu=4*muB; %(T^2*m^3)/J
Ms=mu/V; %T^2/J
Omegae=(2*gamma*abs(4*(-396*kB/V)-532*kB/V))/(mu0*Ms); %s^-1
K4parallel=1.8548*(10^(-25))/V; %J/m^3
Omega4parallel=(2*gamma*K4parallel)/(mu0*Ms); %s^-1
OmegaR=sqrt(2*Omegae*Omega4parallel);
%% Now, let's create the array of times on which we are interested
time=[0:0.001:40].*(10^(-12)); %s
%% Let's define the temperature array
temperature=[0:100];
%% Let's solve the second-order differential equation, where the X variable will have to columns
% Probably, the first column corresponds to phi=x(1) and the second one to
% \dot{phi}=x(2)
lx=zeros(length(time),length(temperature));
ly=zeros(length(time),length(temperature));
mz=zeros(length(time),length(temperature));
for i=1:length(temperature)
    [T,X]=ode45(@(t,x)fx(t,x,Hx,Hy,OmegaR(i),Omegae(i),gamma(i),alpha(i)),time,initial_conditions);
    lx(:,i)=cos(X(:,1));
    ly(:,i)=sin(X(:,1));
    mz(:,i)=-2*Omegae*X(:,2);
    u1=figure(1)
    plot(T.*(10^(12)),lx(:,i),'LineWidth',2)
    hold on
    xlabel('$t \, \, \left( \mathrm{ps} \right)$','FontSize',12,'interpreter','latex')
    ylabel('$l_{\mathrm{x}}$','FontSize',12,'interpreter','latex')
    l1=legend('$T=$',num2str(temperature(i)),' K','FontSize',12,'Location','best');
    set(l1,'interpreter','latex','FontSize',16)
    set(gca,'TickLabelInterpreter','latex','FontSize',16)
    set(u1,'Units','Inches');
    posu1=get(u1,'Position');
    set(u1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[posu1(3),posu1(4)])
    saveas(gcf,fullfile(outputdir,['lx_Time_Evolution_Depending_on_Temperature.pdf']))
end

Obviously I have not loaded the values of the parameters for each temperature yet, but I want to have the scheme complete before loading, so imagine that I have the values already. I have problems on how to implement the legend properly and my approach will save the final plot or not (I would said that my approach won't save the final result because it is inside the for loop).

[Star Strider]

I am a bit lost.

I would do the integrations in the first loop:

for i=1:length(temperature)
    [T{i},X{i}]=ode45(@(t,x)fx(t,x,Hx,Hy,OmegaR(i),Omegae(i),gamma(i),alpha(i)),time,initial_conditions);
end

so you have them stored in your workspace, then do all the other calculations and plots in a separate, subsequent loop. Also, you only appear to be plotting ‘T’ and ‘lx’. Since ‘time’ is a (1x40001) vector:

time=[0:0.001:40]*1E-12;

(that it would likely be best to shorten), ‘X’ will be a (40001x2) matrix. You can probably eliminate the plotting loop entirely by plotting the ‘X’ values as matrices, because they will all be the same lengths. You can then do something like this:

Xmat = cell2mat(X);
lx = cos(Xmat(:,1:2:end));
ly = sin(Xmat(:,1:2:end));
mz = -2*Omegae*Xmat(:,2:2:end);

u1=figure(1)
plot(T*1E12,lx,'LineWidth',2)
hold on
xlabel('$t \, \, \left( \mathrm{ps} \right)$','FontSize',12,'interpreter','latex')
ylabel('$l_{\mathrm{x}}$','FontSize',12,'interpreter','latex')
l1=legend('$T=$',num2str(temperature(i)),' K','FontSize',12,'Location','best');
set(l1,'interpreter','latex','FontSize',16)
set(gca,'TickLabelInterpreter','latex','FontSize',16)
set(u1,'Units','Inches');
posu1=get(u1,'Position');
set(u1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[posu1(3),posu1(4)])
saveas(gcf,fullfile(outputdir,['lx_Time_Evolution_Depending_on_Temperature.pdf']))

I tested that as best I could with some random matrices, and it works. You need to use it with your code to test it specifically.

[Richard Wood]

Hey! Thank you very much. I will test it soon in the following days. For the moment, I have only written the plot of lx(t), but it will also implement the one concerning ly(t) and mz(t), but they are completely equivalent to the one for lx(t).

Nevertheless, I have some conceptual questions about your last response:

• The element X{i}, as you have said, has two columns. You said that, due to the number of elements of time, X must be 40001x2 elements. I am assuming that the first column of X corresponds to x(1) defined in my function file (that it is, to phi), and that the second one corresponds to x(2), that it is, the time derivative of x(1). This is correct? Moreover, after the loop with temperatures, as far as I can imagine, T and X will be structured in the following form: two three-dimensional arrays, being T 40001x1xlength(temperature) and X being 40001x2xlength(temperature). This is true?
• Seeing that after that you use cell2mat, I am more confused. Why do you use that function and which is the aspect of the arrays X and T after the for loop?
• Taking into account this, should I delete what I put previously, that it is, define before enter in the for loop the lx, ly and mz matrices as follows?

lx=zeros(length(time),length(temperature));
ly=zeros(length(time),length(temperature));
mz=zeros(length(time),length(temperature));

• In the same line, I am not sure about the syntaxis present in your previous comment:

lx = cos(Xmat(:,1:2:end));
ly = sin(Xmat(:,1:2:end));
mz = -2*Omegae*Xmat(:,2:2:end);

For sure it is OK, but I am trying to understand it. If I see the structure of index implicit in the use of Xmat, it seems that it is a two-dimensional array. However, I do not understand why you put 1:2:end (probably because I do not understand out which is the precise outcome of X, T. For me, it is like you were taking all the rows of a two-dimensional array, but jumping to the consecutive odd position of all the columns. I suppose that at the end you have as much pairs of zeros(length(time),length(temperature)) matrices as times imposed by the for loop, right?

• Moreover, I need to introduce in the for loop, probably, mz to be calculated, because mz depends on Omegae(i). This should be a problem due to your proposal for calculate lx, ly and mz out of the for loop? I mean, you define out of the loop the Xmat, which is what you use to evaluate mz
• Finally, I am a little bit confused with a part of the plot. The inclusion of a hold on after the syntaxis of plot makes the impression that you are plotting lx(t) for all the temperature possibilities as if you had a for loop. I mean, for each temperature value, a curve lx(t). Moreover, I suppose that, concerning the legend, you have not paid attention to it and that won't work, right? For me to know MatLab what I want to appear num2str(temperature(i)), it should know exactly on which temperature value I am. I can only think about that in terms of a for loop again (sorry, maybe I'm not very modern to properly exploit the potential offered by MatLab)

Thank you very much for your attention.

[Star Strider]

As always, my pleasure!

1) I would delete the current preallocations. I would add one for the cell array:

X = cell(1, length(temperature);

2) The first column of ‘X’ corresponds to ‘x(1)’ in ‘fx’ and the second column to ‘x(2)’.

3) I use cell2mat to create a matrix out of the cell array. Saving to a cell array in the loop is easier than creating a double matrix and concatenating them. The matrix is then composed of:

i    =     1        2     ...
Xmat = [X(:,1:2) X(:,1:2) ... ]

so the ‘X(:,1)’ and ‘X(:,2)’ values alternate column-wise, with all the ‘X(:,1)’ values being ‘Xmat(:,[1 3 5 ...])’, and the ‘X(:,2)’ values being ‘Xmat(:,[2 4 6 ...])’. I doubt that any indexing scheme would be entirely straightforward here.

4) ‘... I need to introduce in the for loop, probably, mz to be calculated, because mz depends on Omegae(i).’

The ‘Omegae’ variable was not originally subscripted. If ‘Omegae’ is a vector the same column-size as ‘X’ (and you have R2016b or later), you can take advantage of ‘automatic implicit expansion’ to calculate ‘mz’ without an explicit loop:

mz = -2*Omegae.*Xmat(:,2:2:end);

If you have an earlier version, you need to use bsxfun for this multiplication:

mz = -2*bsxfun(@times, Omegae, Xmat(:,2:2:end));

The loops are all still there, just much more efficient than explicit for looops.

5) With respect to the plot, I just copied your code and only made changes to the original plot arguments, nothing else. The hold call is not necessary if you are plotting a matrix.

6) With respect to the legend, to make it compatible with my code, it becomes:

l1=legend(compose('$T = %.1f$ K', temperature),'FontSize',12,'Location','best');

Make any necessary changes to it to get the result you want.

I beleive that covers everything. If I omitted anything or if anything needs to be clarified, I will do my best to provide them.

[Richard Wood]

Today I had the opportunity of running the code that you had helped me to implement, being its final form the following:

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];
% Also, I need to define a parameter, Omegae
Neel_Temperature=1575; %K
beta=0.333;
mu0=4*pi*10^(-7); %H/m
gamma=2.21*10^5; %m/(A*s)
kB=1.38064852*10^(-23); %J/K
a0=3.328*10^(-10); %m
c=8.539*10^(-10); %m
V=c*(a0^2); %m^3
Hx=0;
Hy=40*79.5774715459; % Oe to A/m
muB=9.27400994*10^(-24); %(T^2*m^3)/J
mu=4*muB; %(T^2*m^3)/J
Ms=mu/V; %T^2/J
K4parallel=1.8548*(10^(-25))/V; %J/m^3
alpha_T0=0.001;
Omegae_T0=(2*gamma*abs(4*(-396*kB/V)-532*kB/V))/(mu0*Ms); %s^-1
Omega4parallel_T0=(2*gamma*K4parallel)/(mu0*Ms); %s^-1
OmegaR_T0=sqrt(2*Omegae_T0*Omega4parallel_T0);
%% Now, let's calculate the parameters that depend on temperature
Temperature=[0:200:1400]; %K
magnetization_temperature=(1-(Temperature./Neel_Temperature)).^(beta);
alpha=alpha_T0.*(1-(Temperature./(3*Neel_Temperature)));
Omegae=Omegae_T0.*magnetization_temperature; % s^(-1)
OmegaR=OmegaR_T0.*(magnetization_temperature.^5); % s^(-1)
%% Now, let's create the array of times on which we are interested
time=[0:0.001:40].*(10^(-12)); %s
%% Let's solve the second-order differential equation, where the X variable will have to columns
% Probably, the first column corresponds to phi=x(1) and the second one to
% \dot{phi}=x(2)
for i=1:length(Temperature)
    [T{i},X{i}]=ode45(@(t,x)fx(t,x,Hx,Hy,OmegaR(i),Omegae(i),gamma,alpha(i)),time,initial_conditions);
end
Xmat=cell2mat(X);
Tmat=cell2mat(T);
lx=cos(Xmat(:,1:2:end));
ly=sin(Xmat(:,1:2:end));
mz=-2*Omegae.*Xmat(:,2:2:end);
u1=figure(1)
plot(Tmat.*(10^(12)),lx,'LineWidth',2)
xlabel('$t \, \, \left( \mathrm{ps} \right)$','FontSize',14,'interpreter','latex')
ylabel('$l_{\mathrm{x}}$','FontSize',14,'interpreter','latex')
t1=title('Switching process and its dependence on temperature','FontSize',14,'interpreter','latex')
set(t1,'interpreter','latex','FontSize',14)
l1=legend(compose('$T = %.1f$ K',Temperature),'FontSize',14,'Location','northeastoutside');
set(l1,'interpreter','latex','FontSize',14)
set(gca,'TickLabelInterpreter','latex','FontSize',14)
set(u1,'Units','Inches');
posu1=get(u1,'Position');
set(u1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[posu1(3),posu1(4)])
saveas(gcf,fullfile(outputdir,['lx_Time_Evolution_Depending_on_Temperature.pdf']))

Apparently there is no error, because I obtain something without warnings. Nevertheless, it would be great to implement a new thing. That it is, I had the function defined as:

   function dx=fx(t,x,Hx,Hy,OmegaR,Omegae,gamma,alpha);
     % Here, x represents the in-plane angle, phi, and t represents the
     % time variable
     % Global variables: Hx, Hy, OmegaR, Omegae, gamma, alpha
     % Here, x(1)=phi, x(2)=\dot{phi}
     dx=zeros(2,1);
     dx(1)=x(2);
     dx(2)=-(OmegaR.^2/4).*sin(4*x(1))-2*Omegae*gamma*(Hy*cos(x(1))-Hx*sin(x(1)))+...
         2*Omegae*alpha*x(2);
   end

The thing is that I want to implement the same, but depending in this case Hy (defined above as a constant on time and temperature) on time. From time=0 to 5*10^(-12), it will be constant and have the value Hy=40*79.5774715459. For the rest of the times, this is going to be Hy = 0. Do you have some advice and on which of the two scripts I have to focus more to introduce this temporal dependence on Hy?

Obviously it is easy to write the function for Hy on time:

for i=1:length(time)
   if time(i)<=(5.*(10^(-12)))
       Hy(i)=40*79.5774715459; % Oe to A/m
   else Hy(i)=0;
   end
end

Sorry for asking an additional question relation to this.

[Star Strider]

I would set it up as an anonymous function:

Hy_fcn = @(t) (40*79.5774715459).*(t <= 5E-12);

(assuming ‘t’ will never be less than 0), then call it as:

Hy = Hy_fcn(t);

You can call it anywhere in your ‘fx’ function after you define it.

The one caution is that this creates a ‘step’ discontinuity that is not itself differentiable, and will cause problems for any numerical integration routine. It would be best to integrate up to 5E-12, with ‘Hy’ defined as 40*79.5774715459, save the last result of the previous integration as the new initial conditions (and time), then re-start the integration with ‘Hy’ defined as 0, and the new ‘time’ vector beginning at 5E-12. This requires a loop, however it willl give the correct result. Vertically concatenate the ‘T’ and ‘X’ arrays (respectively) to get the full result.

If ‘t’ can be anything, add an extra logical condition to restrict it to being nonzero only in that region:

Hy_fcn = @(t) (40*79.5774715459).*((t >= 0) & (t <= 5E-12));

It will always be 0 outside that region (regardless of which version of the function you use), so an additional term definiing it as such is not necessary.

Note that:

5*10^(-12)

requires a computation, while:

5E-12

does not, and is also easier to read.

[Richard Wood]

Hey! Thank you for your prompt response! First, thanks for pointing out that the way I wrote 5E-12 was computationally slower. Also, thanks for pointing out the discontinuity due to Hy, I had it in mind but I hadn't tried yet to implement it.

From your previous comment, and from what I have understood of it and what I have been thinking later, I have done the following. The script that contains the fx function has remained unchanged. With respect to the other, I have first defined two Hy, Hy1 and Hy2 with the two interest values. Then, I looked for the index that occupies in the array time the last element <= 3E-12 (it was 3E-12, not 5E-12, I was wrong in the previous comment). For some reason that find fails, but the idea is clear. Find the index and then create two temporary arrays that contain the interest values. Once this is done, I have implemented two for loops for integration using ode45 that are equivalent, except for the different time arrays, Hy's and initial conditions. The point is that I am not sure how to write in the elements of the X and T arrays the correct syntax to find the initial condition that I will use for each temperature in the second loop (as you see, I left it blank for now). What do you think?

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_conditions1=[xinitial dxinitial];
% Also, I need to define a parameter, Omegae
Neel_Temperature=1575; %K
beta=0.333;
mu0=4*pi*10^(-7); %H/m
gamma=2.21*10^5; %m/(A*s)
kB=1.38064852*10^(-23); %J/K
a0=3.328*10^(-10); %m
c=8.539*10^(-10); %m
V=c*(a0^2); %m^3
Hx=0;
Hy1=40*79.5774715459;
Hy2=0;
muB=9.27400994*10^(-24); %(T^2*m^3)/J
mu=4*muB; %(T^2*m^3)/J
Ms=mu/V; %T^2/J
K4parallel=1.8548*(10^(-25))/V; %J/m^3
alpha_T0=0.001;
Omegae_T0=(2*gamma*abs(4*(-396*kB/V)-532*kB/V))/(mu0*Ms); %s^-1
Omega4parallel_T0=(2*gamma*K4parallel)/(mu0*Ms); %s^-1
OmegaR_T0=sqrt(2*Omegae_T0*Omega4parallel_T0);
%% Now, let's calculate the parameters that depend on temperature
Temperature=0; %K
magnetization_temperature=(1-(Temperature./Neel_Temperature)).^(beta);
alpha=alpha_T0.*(1-(Temperature./(3*Neel_Temperature)));
Omegae=Omegae_T0.*magnetization_temperature; % s^(-1)
OmegaR=OmegaR_T0.*(magnetization_temperature.^5); % s^(-1)
%% Now, let's create the array of times on which we are interested
time=[0:0.001:40].*(10^(-12)); %s
critical_index_time=find(time<=3E-12,'last');
time_first_interval=time(1:critical_index_time);
time_second_interval=time((critical_index_time+1):end);
%% Let's solve the second-order differential equation, where the X variable will have to columns
% Probably, the first column corresponds to phi=x(1) and the second one to
% \dot{phi}=x(2)
for i=1:length(Temperature)
    [T1{i},X1{i}]=ode45(@(t,x)fx(t,x,Hx,Hy1,OmegaR(i),Omegae(i),gamma,alpha(i)),time_first_interval,initial_conditions1);
end
Xmat1=cell2mat(X1);
Tmat1=cell2mat(T1);
initial_conditions2=[Xmat1() Xmat1()];
for i=1:length(Temperature)
    [T2{i},X2{i}]=ode45(@(t,x)fx(t,x,Hx,Hy2,OmegaR(i),Omegae(i),gamma,alpha(i)),time_second_interval,initial_conditions2);
end
Xmat2=cell2mat(X2);
Tmat2=cell2mat(T2);
Xmat=[Xmat1 Xmat2];
Tmat=[Tmat1 Tmat2];
lx=cos(Xmat(:,1:2:end));
ly=sin(Xmat(:,1:2:end));
mz=-2*Omegae.*Xmat(:,2:2:end);
u1=figure(1)
plot(TmatE-12,lx,'LineWidth',2)
xlabel('$t \, \, \left( \mathrm{ps} \right)$','FontSize',14,'interpreter','latex')
ylabel('$l_{\mathrm{x}}$','FontSize',14,'interpreter','latex')
t1=title('$\alpha=0.001$, H^{\mathrm{SO}}_{\mathrm{y}}=40 \, \mathrm{Oe}, \tau_{\mathrm{p}}=3 \, \mathrm{ps}','FontSize',14,'interpreter','latex')
set(t1,'interpreter','latex','FontSize',14)
l1=legend(compose('$T = %g$ K',Temperature),'FontSize',14,'Location','northeastoutside');
set(l1,'interpreter','latex','FontSize',14)
set(gca,'TickLabelInterpreter','latex','FontSize',14)
set(u1,'Units','Inches');
posu1=get(u1,'Position');
set(u1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[posu1(3),posu1(4)])
saveas(gcf,fullfile(outputdir,['lx_Time_Evolution_Depending_on_Temperature.pdf']))
u2=figure(2)
plot(TmatE-12,ly,'LineWidth',2)
xlabel('$t \, \, \left( \mathrm{ps} \right)$','FontSize',14,'interpreter','latex')
ylabel('$l_{\mathrm{y}}$','FontSize',14,'interpreter','latex')
t2=title('Switching process and its dependence on temperature','FontSize',14,'interpreter','latex')
set(t2,'interpreter','latex','FontSize',14)
l2=legend(compose('$T = %g$ K',Temperature),'FontSize',14,'Location','northeastoutside');
set(l2,'interpreter','latex','FontSize',14)
set(gca,'TickLabelInterpreter','latex','FontSize',14)
set(u2,'Units','Inches');
posu2=get(u2,'Position');
set(u2,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[posu2(3),posu2(4)])
saveas(gcf,fullfile(outputdir,['ly_Time_Evolution_Depending_on_Temperature.pdf']))
u3=figure(3)
plot(TmatE-12,mz,'LineWidth',2)
xlabel('$t \, \, \left( \mathrm{ps} \right)$','FontSize',14,'interpreter','latex')
ylabel('$m_{\mathrm{z}}$','FontSize',14,'interpreter','latex')
t3=title('Switching process and its dependence on temperature','FontSize',14,'interpreter','latex')
set(t3,'interpreter','latex','FontSize',14)
l3=legend(compose('$T = %g$ K',Temperature),'FontSize',14,'Location','northeastoutside');
set(l3,'interpreter','latex','FontSize',14)
set(gca,'TickLabelInterpreter','latex','FontSize',14)
set(u2,'Units','Inches');
posu3=get(u3,'Position');
set(u3,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[posu3(3),posu3(4)])
saveas(gcf,fullfile(outputdir,['mz_Time_Evolution_Depending_on_Temperature.pdf']))

[Star Strider]

This is difficult for me to read.

If you know that you want to integrate your system only up to 3E-12 with one value for ‘Hy’ and beyond that with another value for ‘Hy’, I would break it upt into a two-step loop. With respect to the integration time vectors, I would use:

time1 - linspace(0, 3E-12, 50);
time2 = linspace(3E-12, 4E-11, 150);

or whatever lengths you want them, defining ‘Hy’ (possibly using ‘Hy_fcn’, since it would not change during the integration) for each segment.

Then, vertically concatenate the ‘T’ vectors and the ‘X’ matrices to get the complete result, returning them in their cell arrays, for each value of ‘temperature’.

[Richard Wood]

Thank you again for your comment. With my previous answers I was not trying to search for a difficult route to solve the problem or to avoid your hints, but I have tried to implement your ideas in a way that I can understand. First, the array time will not present any problem related to values lower than zero. You had proposed to create a function on time to describe the evolution of Hy, let's say

Hy_fcn = @(t) (40*79.5774715459).*(t <= 3E-12);

and then call it as

Hy = Hy_fcn(t);

What I understood from here it is that you are creating an array for Hy in the range 0 to 3E-12, where it stays as a constant value that I have proposed in the past. Nevertheless, I do not understand the advantage of this. I mean, you have to deal with two integration regions, for two time arrays, being Hy constant but with a different value. Moreover, I am not sure if at the end I am able to understand your approach where to define this Hy_fcn function. I suppose that in my fx function script, before the of x(1) and x(2).

After that, you propose to define the time two time arrays using linespace instead of [ : ], okay. But then you propose to define Hy (possibly using ‘Hy_fcn’, since it would not change during the integration) for each segment. Sorry, but I have no idea of what are you trying to say here. Of course you are saying that Hy has a different value on each of those time regions, and accordingly I have to deal with its definition somewhere. Of course I agree. But I do not understand how Hy_fcn it is useful here or not even how to deal with it the scheme of the ode45 inside of the for loop.

I suppose also that you do not like at all my approach to try to implement the two for loops to deal with the integration on each region:

for i=1:length(Temperature)
    [T1{i},X1{i}]=ode45(@(t,x)fx(t,x,Hx,Hy1,OmegaR(i),Omegae(i),gamma,alpha(i)),time_first_interval,initial_conditions1);
end
Xmat1=cell2mat(X1);
Tmat1=cell2mat(T1);
initial_conditions2=[Xmat1() Xmat1()];
for i=1:length(Temperature)
    [T2{i},X2{i}]=ode45(@(t,x)fx(t,x,Hx,Hy2,OmegaR(i),Omegae(i),gamma,alpha(i)),time_second_interval,initial_conditions2);
end
Xmat2=cell2mat(X2);
Tmat2=cell2mat(T2);
Xmat=[Xmat1 Xmat2];
Tmat=[Tmat1 Tmat2];

I think that this approach is not very far from what you are proposing, even when I am not sure on how you would implement it. Nevertheless, as I said, if this scheme is not far from yours, I am dealing with the main problem of define the initial conditions for x(1) and x(2), for each temperature, taking the final value of the integration process of the previous for loop. I imagine that if you were who is implementing this, due to the plan that you propose, you will not need to define the initial conditions as me, but in other, more clever, way.

So to summarize. (1) How the Hy_fcn function will work in each of the time regimes if I define it as you propose, as a function of t? (2) Do you think that the plan of separating the integration regimes in two "indepent" for loops is far from your idea? Would your scheme be more correlated, in terms of the connection between the two for loops? (3) On the definition of the initial condition for the region on which Hy is zero, would it be something like...

initial_conditions2=[Xmat1(end,1:2:end) Xmat1(end,2:2:end)]

For sure the definition of this initial conditions should be inside the (first) Temperature loop, and how to deal with each value corresponding to each temperature seems a challenge.

Sorry for the multiple comments and questions again. I would like to support you more on this page, but it seems that I cannot much that liking your comment and accepting it as an answer. If I can do it on other way, let me know.

[Star Strider]

‘Nevertheless, I do not understand the advantage of this.’

I thouight I expl;ained it in this Comment. Quoting:

‘The one caution is that this creates a ‘step’ discontinuity that is not itself differentiable, and will cause problems for any numerical integration routine. It would be best to integrate up to 5E-12, with ‘Hy’ defined as 40*79.5774715459, save the last result of the previous integration as the new initial conditions (and time), then re-start the integration with ‘Hy’ defined as 0, and the new ‘time’ vector beginning at 5E-12. This requires a loop, however it willl give the correct result. Vertically concatenate the ‘T’ and ‘X’ arrays (respectively) to get the full result.’

‘Sorry, but I have no idea of what are you trying to say here.’

I was simply suggesting two ways of defining ‘Hy’ for the two regions of integration. You can do it however you wish.

‘I suppose also that you do not like at all my approach to try to implement the two for loops to deal with the integration on each region:’

I have no opinion on that. I am simply suggesting efficient and mathematically correct approaches to solving your problem. If it works, and is reasonably efficient, the exact method is not relevant. Coding styles differ.

I was proposing a separate for loop of two steps to account for the changes in ‘Hy’ for differing time regimens, since abruptly changing it within the integration will cause problems for the numerical integrators.

[Richard Wood]

I understood perfeclty your approach and concerns, sorry if from my previous message that was not clear. Anyway, I have been able to implement it, based on your guidance but on my own terms. That closes the question entirely. Thank for your useful answers!

[Star Strider]

As always, my pleasure!