ode45 inside a nested loop

Hello everyone,

My situation is the following. I want to solve a second-order differential equation, where the parameters involved depends on external parameters (in my case, on variables Temperature and Hy1_value). This parameters are Omegae, OmegaR, alpha. My function file looks as follows

   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

Now, I want to solve my second-order differential equation which is inside a nested for loop. The thing is that I am not sure of how to index the variables involved in the second-order differential equation solving process by using ode45. My code looks as follows:

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;
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);
pulse_duration=1E-12; % s
decay_time=100E-12; % s
time=[0:0.001:70].*(10^(-12)); % s
critical_index_time=find(time==1E-12);
time_first_interval=time(1:critical_index_time); % s
time_second_interval=time((critical_index_time+1):end); % s
Temperature_max=[200:200:1200]; % K
Hy1_values=linspace(100,1000,10).*79.5774715459; % A/m
Temperature1=zeros(length(time_first_interval),length(Temperature_max));
Temperature2=zeros(length(time_second_interval),length(Temperature_max));
Temperature=zeros(length(time),length(Temperature_max));
Hy1=zeros(length(time_first_interval),length(Hy1_values),length(Temperature_max));
magnetization_temperature=zeros(length(time),length(Temperature_max));
magnetization_temperature1=zeros(length(time_first_interval),length(Temperature_max));
magnetization_temperature2=zeros(length(time_second_interval),length(Temperature_max));
alpha=zeros(length(time),length(Temperature_max));
Omegae=zeros(length(time),length(Temperature_max)); % s^(-1)
OmegaR=zeros(length(time),length(Temperature_max)); % s^(-1)
alpha1=zeros(length(time_first_interval),length(Temperature_max));
Omegae1=zeros(length(time_first_interval),length(Temperature_max)); % s^(-1)
OmegaR1=zeros(length(time_first_interval),length(Temperature_max)); % s^(-1)
alpha2=zeros(length(time_second_interval),length(Temperature_max));
Omegae2=zeros(length(time_second_interval),length(Temperature_max)); % s^(-1)
OmegaR2=zeros(length(time_second_interval),length(Temperature_max)); % s^(-1)
Hy1=zeros(length(time_first_interval),length(Hy1_values),length(Temperature_max)); % A/m
Hy2=zeros(length(time_second_interval),length(Hy1_values),length(Temperature_max));
Xmat1=zeros(length(time_first_interval),2*length(Hy1_values),length(Temperature_max));
%% Now, let's calculate the parameters that depend on temperature
for i=1:length(Temperature_max)
    for j=1:length(Hy1_values)
        Temperature1(:,i)=(Temperature_max(i)./sqrt(pulse_duration)).*sqrt(time_first_interval); % K
        Temperature2(:,i)=Temperature1(end,i).*(exp(-time_second_interval./pulse_duration)); % K
        Temperature(:,i)=vertcat(Temperature1(:,i),Temperature2(:,i));% K
        magnetization_temperature1(:,i)=(1-(Temperature1(:,i)./Neel_Temperature)).^(beta);
        magnetization_temperature2(:,i)=(1-(Temperature2(:,i)./Neel_Temperature)).^(beta);
        magnetization_temperature(:,i)=vertcat(magnetization_temperature1(:,i),...
            magnetization_temperature2(:,i));
        alpha1(:,i)=alpha_T0.*(1-(Temperature1(:,i)./(3*Neel_Temperature)));
        alpha2(:,i)=alpha_T0.*(1-(Temperature2(:,i)./(3*Neel_Temperature)));
        alpha(:,i)=vertcat(alpha1(:,i),alpha2(:,i));
        Omegae1(:,i)=Omegae_T0.*magnetization_temperature1(:,i); % s^(-1)
        Omegae2(:,i)=Omegae_T0.*magnetization_temperature2(:,i); % s^(-1)
        Omegae(:,i)=vertcat(Omegae1(:,i),Omegae2(:,i)); % s^(-1)
        OmegaR1(:,i)=OmegaR_T0.*(magnetization_temperature1(:,i).^5); % s^(-1)
        OmegaR2(:,i)=OmegaR_T0.*(magnetization_temperature2(:,i).^5); % s^(-1)
        OmegaR(:,i)=vertcat(OmegaR1(:,i),OmegaR2(:,i));
        Hy1(:,j,i)=ones(length(time_first_interval),1).*Hy1_values(j); % A/m
        Hy2(:,j,i)=zeros(length(time_second_interval),1); % A/m
        Hy(:,j,i)=vertcat(Hy1(:,j,i),Hy2(:,j,i));
        %% 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)
        % I DON'T KNOW HOW TO DEAL WITH THE INDEXING FROM HERE ON!
        [T1{i,j},X1{i,j}]=ode45(@(t,x)fx(t,x,Hx,Hy1(:,j,i),OmegaR1(:,i),Omegae1(:,i),gamma,alpha1(:,i)),time_first_interval,initial_conditions1);
%         Xmat1(:,2*j-1,i)=cell2mat(X1);
%         Tmat1=cell2mat(T1);
%         [T2,X2]=ode45(@(t,x)fx(t,x,Hx,Hy2,OmegaR,Omegae,gamma,alpha),time_second_interval,[Xmat1(end,1) Xmat1(end,2)]);
%         Xmat2=cell2mat(X2);
%         Tmat2=cell2mat(T2);
%         Xmat=vertcat(Xmat1,Xmat2);
%         Tmat=vertcat(Tmat1,Tmat2);
%         lx=cos(Xmat(:,1));
%         ly=sin(Xmat(:,1));
%         mz=-(1./(2*Omegae)).*Xmat(:,2);
    end
end

I had to separate the solving process in two ode45's because of the abrupt discontinuity imposed in the Hy variable in the intersection of Hy1 and Hy2. The thing is that all variables are well defined. The problem begins when I put the warning "% I DON'T KNOW HOW TO DEAL WITH THE INDEXING FROM HERE ON!". My objetive is that, at the end, my arrays lx, ly and mz have the following dimensions:

lx=zeros(length(time),2*length(Hy1_values),length(Temperature_max));
ly=zeros(length(time),2*length(Hy1_values),length(Temperature_max));
mz=zeros(length(time),2*length(Hy1_values),length(Temperature_max));

So that I have all the information after the for loops well-organized and saved.

Any suggestions?

The error that MatLab gives to me is

Unable to perform assignment because the left and right sides have a different number of elements.

Error in fx (line 8)
     dx(2)=-(OmegaR.^2/4).*sin(4*x(1))+2*Omegae.*gamma.*(Hy.*cos(x(1))-Hx.*sin(x(1)))-...

Error in
Second_Order_Differential_Equation_Laser_Like>@(t,x)fx(t,x,Hx,Hy1(:,j,i),OmegaR1(:,i),Omegae1(:,i),gamma,alpha1(:,i))
(line 80)
        [T1{i,j},X1{i,j}]=ode45(@(t,x)fx(t,x,Hx,Hy1(:,j,i),OmegaR1(:,i),Omegae1(:,i),...

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Second_Order_Differential_Equation_Laser_Like (line 80)
        [T1{i,j},X1{i,j}]=ode45(@(t,x)fx(t,x,Hx,Hy1(:,j,i),OmegaR1(:,i),Omegae1(:,i),...

The thing is that OmegaR, Omegae and alpha depends on time. I have seen that when you are dealing time-dependence of the parameters involved on your differential equation, when using ode45, you must put (t) after then to reflect its dependence on time, but I think that I should also take care of the indexing of that multidimensional arrays!

[darova]

Here is the problem

I suppose there is some dependence between Hx and t (or x)?

Can you show original equations/formulas,

I also suggest you to give shorter names to your variables. Two weeks will pass and you will not be able to read your code

t1int = time(1:critical_index_time); % s
t2int = time((critical_index_time+1):end); % s
Temp_max=[200:200:1200]; % K
[Temp1, mag_temp1, alpha1, Omegae1, OmegaR1] = deal( zeros(length(t1int),length(Temp_max)) );
[Temp2, mag_temp2, alpha2, Omegae2, OmegaR2] = deal( zeros(length(t2int),length(Temp_max)) );
[Temp, mag_temp, alpha, Omegae, OmegaR] = deal( zeros(length(time),length(Temp_max)) );

[Roderick]

Hey! Thank for your answer. I have simplified a little bit one of the scripts, as you have proposed. I hope that that helps to read it easily. Again, my function script:

   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 other script:

%% 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 certain parameters of my material
NeelT=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;
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);
pulse_duration=1E-12; % s
decay_time=100E-12; % s
%% Now,I can define the time array and also the parameters that will depend
% on Temperature (Temp, T1int, T2int). That it is, OmegaR(1,2),
% Omegae(1,2), alpha(1,2). The thing is that temperature will depend on
% time, so these coefficients that depend on temperature will depend on
% time
time=[0:0.001:70].*(10^(-12)); % s
critical_index_time=find(time==1E-12);
t1int=time(1:critical_index_time)'; % s
t2int=time((critical_index_time+1):end)'; % s
Tmax=[200:200:1200]; % K
Hy1_values=linspace(100,1000,10).*79.5774715459; % A/m
Hy1=zeros(length(t1int),length(Hy1_values),length(Tmax)); % A/m
Hy2=zeros(length(t2int),length(Hy1_values),length(Tmax));
Hy=zeros(length(time),length(Hy1_values),length(Tmax));
[magT1, Omegae1, OmegaR1, alpha1, T1int]=deal(zeros(length(t1int),length(Tmax)));
[magT2, Omegae2, OmegaR2, alpha2, T2int]=deal(zeros(length(t2int),length(Tmax)));
[magT, Omegae, OmegaR, alpha, Temp]=deal(zeros(length(time),length(Tmax)));
%% Now, let's calculate the parameters that depend on temperature (indirectly on time)
for i=1:length(Tmax)
    for j=1:length(Hy1_values)
        T1int(:,i)=(Tmax(i)./sqrt(pulse_duration)).*sqrt(t1int);
        T2int(:,i)=T1int(end,i).*(exp(-t2int./pulse_duration));
        Temp(:,i)=vertcat(T1int(:,i),T2int(:,i));
        magT1(:,i)=(1-(T1int(:,i)./NeelT)).^(beta);
        magT2(:,i)=(1-(T2int(:,i)./NeelT)).^(beta);
        magT(:,i)=vertcat(magT1(:,i),magT2(:,i));
        alpha1(:,i)=alpha_T0.*(1-(T1int(:,i)./(3*NeelT)));
        alpha2(:,i)=alpha_T0.*(1-(T2int(:,i)./(3*NeelT)));
        alpha(:,i)=vertcat(alpha1(:,i),alpha2(:,i));
        Omegae1(:,i)=Omegae_T0.*magT1(:,i);
        Omegae2(:,i)=Omegae_T0.*magT2(:,i);
        Omegae(:,i)=vertcat(Omegae1(:,i),Omegae2(:,i)); % s^(-1)
        OmegaR1(:,i)=OmegaR_T0.*(magT1(:,i).^5); % s^(-1)
        OmegaR2(:,i)=OmegaR_T0.*(magT2(:,i).^5); % s^(-1)
        OmegaR(:,i)=vertcat(OmegaR1(:,i),OmegaR2(:,i));
        Hy1(:,j,i)=ones(length(t1int),1).*Hy1_values(j); % A/m
        Hy2(:,j,i)=zeros(length(t2int),1); % A/m
        Hy(:,j,i)=cat(2,[Hy1(:,j,i);Hy2(:,j,i)]);
        %% Let's solve the second-order differential equation, where the X variable will have to columns
        [T1{j},X1{j}]=ode45(@(t,x)fx(t,x,Hx,Hy1,OmegaR1(:,i),Omegae1(:,i),...
            gamma,alpha1(:,i)),t1int,initial_conditions1);
% %         Xmat1(:,2*j-1,i)=cell2mat(X1);
% %         Tmat1=cell2mat(T1);
% %         [T2,X2]=ode45(@(t,x)fx(t,x,Hx,Hy2,OmegaR2,Omegae2,gamma,alpha2),t2int,[Xmat1(end,1) Xmat1(end,2)]);
% %         Xmat2=cell2mat(X2);
% %         Tmat2=cell2mat(T2);
% %         Xmat=vertcat(Xmat1,Xmat2);
% %         Tmat=vertcat(Tmat1,Tmat2);
% %         lx=cos(Xmat(:,1));
% %         ly=sin(Xmat(:,1));
% %         mz=-(1./(2*Omegae)).*Xmat(:,2);
    end
end

I am sorry, but I do not see the point on your previous kind comment. Which is exactly the problem with the indexing? In my problem Hx=0, always. I suppose that your previous questions refers to Hy. The point is the following. I want to solve that second-order differential equation for different temperature profiles (external for loop) for a variety of magnetic fields (internal for loop, Hy variable, which values are enconded in Hy1_values, because in the frontier between t1int and t2int Hy goes from a constant value to be zero abruptly, that is the reason of using two ode45 solvers between these two problematic regions). Hy depends on time, of course. I am not sure that I need to define such a difficult three-dimensional array for Hy. I only want to play with different constant values for Hy on the t1int for different temperature profiles.

At the end what I want to do is to plot lx, ly and mz. One plot for each value of Tmax, being on each of them all the possible cases for Hy1.

I attach a file with the formulas in LaTeX (with not a great format, but it is better than look into MatLab syntaxis). I hope this helps a little bit!

Answer 1

[darova]

I want to show you a simple example of how ode45 works. I wrote simple solver like ode45

function main
clc
F = @(t,x) [x(2); x(1)-2*sin(x(1))];
ts = [0 5];
x0 = [0.5 1];
[t1,x] = ode45(F,ts,x0);                    % standard ode45 solver
[t2, res] = new_ode45_solver(F,ts,x0);      % solve equations with new solver
h1 = plot(t1,x,'r');
hold on
h2 = plot(t2,res,'b');
hold off
legend([h1(1) h2(2)],'ode45 solver','my new ode45 solver')
end

function [t, res] = new_ode45_solver(F,ts,x0)
n = 100;                            % number of steps
dt = (ts(2)-ts(1))/n;               % time step
res = zeros(n,2);                   
t = linspace(ts(1),ts(2),n)';
x = x0;
for i = 1:n
    res(i,:) = x;
    x = x + dt*F(t(i),x)';          % standard Euler scheme/method
end
end

Look at this line

x = x + dt*F(t(i),x)';          % standard Euler scheme/method

Look here

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);
% And look this
x(2) = [1 2 3]              % wring! - doesn't work

ode45 works as loop. So if you pass some parameters (Hx,Hy) they can't be arrays

[darova]

Here is how i see solution for your eqution. Maybe you will be interested

Bigger Hy - red, small Hy - blue

[Roderick]

Hey! What an incredible answer you have given me. I cannot say that I understand everything completely since I am not an expert in MatLab, but I am trying to learn as much as I can.

Regarding your first answer where you have exposed how ode45 really works. First, a silly question. I am surprised that by defining your time array as [0 5], specifying only two values, ode45 is able to understand that you want to solve your differential equation between those two points at a number of intermediate points that it chooses. I didn't know it worked like this, if I understood correctly. Because for example in your case you define your time extremes and define what the time step will be to solve the differential equation. Is the discrepancy between ode45 and your approach somewhat based on the discrepancy between your grids?

Now, regarding you Untitled3.mat file. First of all, thank you very much for getting involved in this with such enthusiasm. I think that maybe I have not explained it very well, or maybe I do not have understood that in your code. A minor details is that you use thau_d for the two parts of the Temp function, when the second is governed by thau_p, not defined in your code. But it does not matter. The thing is Hy0. What I wanted was to made Hy=cte up to time=1E-12, and from that point on time, made it 0 for the rest of solving process. I think that that it is not contemplated in your approach, right? Just to know. Edit: it seems that you have taken into account even that, I have already seen where you have done it, I think, the second function! That syntaxis of saying that Hy is a certain value up to 1E-12 is enough?

I have taken a look and tried to understand what you have done. For me, it is surprising that you structured the problem as a nested function, I had never seen anything like it. I guess if I deleted the main function of the header and the end of the file it would work just as well, right? I have tried and see that it is not so, that it does not work. I will try to see why the need to implement it like this. On the other hand, I am surprised again how to define the vector of times ts, but I see that it works, something new that I learn. What does surprise me is the way you have chosen to define Temp in both time regimes. I imagine that the fact that you put the sum simply implies that when you pass the threshold of 1E-12, the first term of the sum is 0. A very elegant way to deal with the definition of the temperature function. When you enter the two for loop, you define a variable counting, k. However, I do not see that it comes into play anywhere, unless it is being renewed in each cycle. What does the introduction of this variable imply? On the other hand, it is great that you have managed to find a way to solve the differential equation by dealing with only one index in each variable, that is, in Hy0 and Tmax0, that greatly simplifies everything and removes the problem of having to deal with many indices. Regarding your plot on Untitled3.mat, it's great. However, it is not exactly what I was looking for, but thanks for taking the trouble to implement it. What I'm looking for plotting is, for each Tmax0, all cases of Hy0 in a two-dimensional plot. Each case of Tmax0 in a single plot. And for each Tmax0 a plot for lx, another for ly, and another for mz vs time. Also, capture in the title of each plot the value of Tmax0 that is being embodied in that particular plot. Also, let the legend be constructed properly by capturing the Hy0 you are using. Also to save each of these plots. In the past before you implemented this new code, I had already raised a bit how I planned to start doing it:

%         u1=figure(i)
%         plot(Tmat.*(10^(12)),lx,'LineWidth',2)
%         hold on
%         xlabel('$t \, \, \left( \mathrm{ps} \right)$','FontSize',14,'interpreter','latex')
%         ylabel('$l_{\mathrm{x}}$','FontSize',14,'interpreter','latex')
%         xline(1,'--k')
%         xlim([0 time(end).*(10^(12))])
%         t1=title('$\alpha=0.001, \tau_{\mathrm{p}}=1 \, \mathrm{ps}, T_{\mathrm{max}}=$',num2str(Temperature_max(i)),' K','FontSize',14,'interpreter','latex')
%         set(t1,'interpreter','latex','FontSize',14)
%         l1=legend(compose('$H^{\mathrm{SO}}_{y} = %g$ mT',Hy1(1).*(0.1/79.5774715459)),'FontSize',14,'Location','northeast');
%         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_Laser_Like_Tmax_',num2str(Temperature_max(i)),'_K.pdf']))
%         u2=figure(i)
%         plot(Tmat.*(10^(12)),ly,'LineWidth',2)
%         hold on
%         xlabel('$t \, \, \left( \mathrm{ps} \right)$','FontSize',14,'interpreter','latex')
%         ylabel('$l_{\mathrm{y}}$','FontSize',14,'interpreter','latex')
%         xline(1,'--k')
%         xlim([0 time(end).*(10^(12))])
%         t2=title('$\alpha=0.001, \tau_{\mathrm{p}}=1 \, \mathrm{ps}, T_{\mathrm{max}}=$',num2str(Temperature_max(i)),' K','FontSize',14,'interpreter','latex')
%         set(t2,'interpreter','latex','FontSize',14)
%         l2=legend(compose('$H^{\mathrm{SO}}_{y} = %g$ mT',Hy1(1).*(0.1/79.5774715459)),'FontSize',14,'Location','southeast');
%         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_Laser_Like_Tmax_',num2str(Temperature_max(i)),'_K.pdf']))
%         u3=figure(i)
%         plot(Tmat.*(10^(12)),mz,'LineWidth',2)
%         hold on
%         xlabel('$t \, \, \left( \mathrm{ps} \right)$','FontSize',14,'interpreter','latex')
%         ylabel('$m_{\mathrm{z}}$','FontSize',14,'interpreter','latex')
%         xline(1,'--k')
%         xlim([0 time(end).*(10^(12))])
%         t3=title('$\alpha=0.001, \tau_{\mathrm{p}}=1 \, \mathrm{ps}, T_{\mathrm{max}}=$',num2str(Temperature_max(i)),' K','FontSize',14,'interpreter','latex')
%         set(t3,'interpreter','latex','FontSize',14)
%         l3=legend(compose('$H^{\mathrm{SO}}_{y} = %g$ mT',Hy1(1).*(0.1/79.5774715459)),'FontSize',14,'Location','southeast');
%         set(l3,'interpreter','latex','FontSize',14)
%         set(gca,'TickLabelInterpreter','latex','FontSize',14)
%         set(u3,'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_Laser_Like_Tmax_',num2str(Temperature_max(i)),'_K.pdf']))
    end
%     u4=figure(100)
%     plot(time.*(10^(12)),Temperature(:,i),'LineWidth',2)
%     xlabel('$t \, \, \left( \mathrm{ps} \right)$','FontSize',14,'interpreter','latex')
%     ylabel('$T \, \, \left( \mathrm{K} \right)$','FontSize',14,'interpreter','latex')
%     xlim([0 time(end).*(10^(12))])
%     t4=title('Temperature profile for a laser-like excitation','FontSize',14,'interpreter','latex')
%     set(t4,'interpreter','latex','FontSize',14)
%     l4=legend(compose('$T = %g$ K',Temperature_max(i)),'FontSize',14,'Location','best');
%     set(l4,'interpreter','latex','FontSize',14)
%     set(gca,'TickLabelInterpreter','latex','FontSize',14)
%     set(u4,'Units','Inches');
%     posu4=get(u4,'Position');
%     set(u4,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[posu4(3),posu4(4)])
%     saveas(gcf,fullfile(outputdir,['Temperature_Evolution_Laser_Like_Excitation.pdf']))
end

Just for curiosity I was to try to see all the temperature profiles in a same graph, but that it is just a curiosity. I have to think how to implement my idea of the aforementioned plots on your approach. And for that, I had asked you the previous question, to try to understand better the file.

[darova]

• Is the discrepancy between ode45 and your approach somewhat based on the discrepancy between your grids?

It doesn't matter if you use time interval [0 time] or linspace(0,time). ode45 solver choose it's own timestep and then create time span (the same span if you specified). Look HELP about ode45

• The thing is Hy0. What I wanted was to made Hy=cte up to time=1E-12, and from that point on time, made it 0 for the rest of solving process. I think that that it is not contemplated in your approach, right? Just to know.

Sorry. FIrst uploaded file without this. Just add into the function

Hy = (t<1E-12)*Hy;

• What does the introduction of this variable imply?

I just wanted to color each curve in different color

k = 0;
cmap = jet( length(Hy0)*length(Tmax0) );

But later decided to color only by Hy0 changes. So you can remove it

[Roderick]

Thank you for your comments, Darova. I have accepted your answer. If you have any useful comment about the plots, let me know, if you want. I just have some doubt about how LX, LY, MZ values are stored in each loop to undertake my plots. That it is, if I should somehow implement the plots inside the for loop for Hy or if all the values for lx, ly, mz are storaged in that variable after that loop and I can plot them outside for loop of Hy but inside of the one of temperature, of course.