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

Hi.I have a code error.

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

Error in rate_eq_program_1 (line 26)

ys(1) = I./(E.*Vp)-(z(1)./tc)-GN.*(z(1)-Ntr).*sigmas.*((abs(z(3)).^2)./Vc);

%%%%%%%%%%%%%%%%%%%%%%%% constants %%%%%%%%%%%%%%%%%%%%%%%%%
q = 1.6021766208e-19;                       
c = 2.99792458e8;                                                         
kk = 1.38e-23;                                                          
kk_eV = kk/q;                                                            
T = 300;                                                                  
h = 6.626068e-34;                                                         
h_eV = h/q;                                                               
hbar = h/(2*pi);                        
hbar_eV = hbar/q;      
epsilon0 = 8.854187817e-12;  
P = -1; 
rL = 1;
rB = 0.2 ;
tB = sqrt(1-rB^2);
Qt =500;
% Qc =180;
Qi = 14300;
Qp = 10000;
%Qc = 1800;
L = 5e-6;
ng = 3.5;
sigma = 3;
m = 25;
A = 0.105e-12; 
confine = 0.5;
Vp = 0.576e-18; 
Vc = confine*A*L;
gN =5e-20;
Ntr = 1e24;
tc = 0.5e-9; 
%tc = 0.7e-9; 
alphai = 1000;
%alpha = 1.5;
alpha =1;
confinec = 0.3;
vg = c/ng;
vgp = c/6/ng;
Vpc = 0.22e-18;
KD = confinec*vgp*alpha*gN/2; 
KA = -confinec*vgp*gN/2;

lambda0 = 1554e-9;
deltalambda = 0.01e-9;
lambda_start = 1540e-9;
lambda_end = 1560e-9;
lambda = lambda_start:deltalambda:lambda_end;
omega = 2*pi*c./lambda;
omega0 = omega(find(lambda==lambda0));
E = hbar_eV.*omega;
gamai = omega./2/Qi;
gamap = omega./2/Qp;
gamat = omega./2/Qt;
%gamac = omega./2/Qc;
%gamac = omega./2/Qc;
gamac = gamat-gamai-gamap;
gama1 = gamac./2;
gama2 = gama1;
COS2 = 0.5*gama1./gama1.*tB^2/rB-0.5*tB^2/rB-rB;
SIN2 = P*tB./(2.*gama1*rB).*sqrt(4.*gama1.*gama1 - tB^2.*(gamac.^2));
% CS12 = 1/(1i*tB).*(1i*SIN2 + rB);
GLA = gN*vgp./(E.*Vpc*q);
GN = confine*vg*gN;
tin = 2*L/vg;
M = 0.5*(length(E)-1);
Tmin = 0;
Tmax = 5e-9;
dt = 0.01e-12;
Time = Tmin:dt:Tmax;
N0 = zeros(1,4);
N0(3) = 0.2;
%I =2e-3;TA I =12e-3;
I =0.5e-3;
Vnc= 0.24e-18;

%%%%%%%%%%%%%%%%%%%%%%% main code  %%%%%%%%%%%%%%%%%%%%%%%%%%%
Time = Tmin:dt:Tmax;      
z0 = [0.2,0.2,0.2,0.2]; 
[t,z] = ode45(@rate_eq_program_1,Time,z0);
param_rate_eq

figure(1)
plot(t,z(:,1));   % divided to normalize
xlabel('time [ns]','FontSize',14);  % size of x label
ylabel('Arbitrary units','FontSize',14);    % size of y label
set(gca,'FontSize',14);                 % size of tick marks on both axis
legend('N', 'Location','SE')  % legend inside the plot

figure(2)
plot(t,z(:,2));   % divided to normalize
xlabel('time [ns]','FontSize',14);  % size of x label
ylabel('Arbitrary units','FontSize',14);    % size of y label
set(gca,'FontSize',14);                 % size of tick marks on both axis
legend('NC', 'Location','SE')  % legend inside the plot

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function %%%%%%%%%%%%%%%%%%%%%%%%

function ys = rate_eq_program_1(t,z)
param_rate_eq_fano
ys=zeros(4,1);

%rR=(-p*gamma_c)/(1i*(delta_c)+gamma_T-1/2*(1-1i*henry)*conf_NC*V_g*g_n*(z(2)-N_0));
deltawc=KD.*(z(2))-1i*KA.*(z(2));
deltac=omega0+deltawc-omega;
rR=(-P.*gamac)./(1i*(deltac)+gamat-(1/2*(1-1i)*confinec*vg*gN).*(z(2)-Ntr));
rRS=-P.*gamac./(1i.*(deltac)+gamat);
%rRS=rB+(-1i*tB-rB).*gamac./(1i.*(deltac)+gamat);
g=alphai+(1/(2*L)).*log(1./(abs(rRS)).^2);

A=(2*epsilon0*n*ng)./(hbar.*omega0);
B=(1+(abs(rRS))).*(1-(abs(rRS)));
D=((g)-alphai);
H=c./(omega0.*n);
F=imag(rRS)./abs(rRS);
sigmas=A.*((B./D)+(H.*F));


ys(1) = I./(E.*Vp)-(z(1)./tc)-GN.*(z(1)-Ntr).*sigmas.*((abs(z(3)).^2)./Vc);
ys(2) = (-z(2)./tc)-GLA.*(z(2)-Ntr).*(abs(z(4)).^2)./Vnc;
ys(3) = 1/2*(1-1i)*GN.*(z(1)-Nss).*z(3)+(1/tin.*((z(4)./rR)-z(3)));
ys(4) =((-1i.*deltawc-gamat).*z(4))-(P.*gamac.*z(3))+(1/2*(1-1i).*confinec*vg*gN.*(z(2)-Ntr).*z(4));
ys = ys';
end

i tried ys=zeros(4,M) But it's also a mistake.

i tried ys=zeros(4,length(E)) But it's also a mistake.

....

Thanks for helping someone.

[Walter Roberson]

omega0 = omega(find(lambda==lambda0));

You used a floating point equality test on values that are calculated. The == operator does bit-for-bit exact comparison without accepting even a little rounding. Therefore the test will usually come out with no matches.

You should only ever use == with floating point numbers to compare values that have been extracted from the same calculation, such as x==min(x) where you know that a value is bit-for-bit identical. (There are also some good reasons to check for exact 0, and sometimes to compare to integers but less often than you would expect).

Otherwise you should only compare within a tolerance, perhaps using ismembertol()

[mohammad heydari]

Thank you dear Walter.

I understood what you were saying but didn't understand what to do to fix this problem.Can you say exactly the solution to this problem because I have no idea about it.

Best regards

[Ioannis Andreou]

Instead of

lambda=lambda0

use

abs(lambda-lambda0) < 1e-15

[mohammad heydari]

hi.

Problem not resolved

[Walter Roberson]

You should format your code as a block to make it easier for people to copy it.

[Walter Roberson]

What is param_rate_eq_fano ? Is it a script that defines the n used in

    A=(2*epsilon0*n*ng)./(hbar.*omega0);

Answer 1

[Walter Roberson]

Your function rate_eq_program_1 uses a number of variables defined in the main program, but without importing them. You need to add a "function" statement at the very top of your code, and add an "end" statement at the end of your code, in order to make rate_eq_program_1 into a nested function that can access the variables.

Now, in the main part of your program you have

lambda = lambda_start:deltalambda:lambda_end;

so lambda is a vector

omega = 2*pi*c./lambda;

so omega is a vector because it is built from the vector lambda

omega0 = omega(find(lambda==lambda0));

As mentioned before due to floating point round-off the == might not find any matches. It will find either 0 or 1 match, so omega0 will end up either empty or a scalar.

E = hbar_eV.*omega;

E is constructed from the vector omega, so E is a vector. Did you want to construct E from omega0 ?

Now, inside rate_eq_program_1 you have

ys(1) = I./(E.*Vp)-(z(1)./tc)-GN.*(z(1)-Ntr).*sigmas.*((abs(z(3)).^2)./Vc);

As noted above, E is a vector, so the right hand side of the assignment is a vector, but the left hand side only has room for a scalar.

Now look a little further and it turns out that sigmas is a vector as well:

deltac=omega0+deltawc-omega;

You use all of the vector omega, so deltac will be a vector.

rR=(-P.*gamac)./(1i*(deltac)+gamat-(1/2*(1-1i)*confinec*vg*gN).*(z(2)-Ntr));
rRS=-P.*gamac./(1i.*(deltac)+gamat);

Those use the vector deltac so they will both be vectors.

    g=alphai+(1/(2*L)).*log(1./(abs(rRS)).^2);

Uses the vector rRS so g is a vector

    A=(2*epsilon0*n*ng)./(hbar.*omega0);

Use the undefined variable n . If we assume that n is a scalar, then A will be a scalar.

    B=(1+(abs(rRS))).*(1-(abs(rRS)));

Uses the vector rRS so B will be a vector

    D=((g)-alphai);

g is a vector so D will be a vector.

    H=c./(omega0.*n);

Use the undefined variable n. If we assume that n is a scalar, then H will be a scalar.

    F=imag(rRS)./abs(rRS);

Uses the vector rRS so F will be a vector.

    sigmas=A.*((B./D)+(H.*F));

Uses the vector B, the vector D, the vector F, so sigmas will be a vector.

ys(1) = I./(E.*Vp)-(z(1)./tc)-GN.*(z(1)-Ntr).*sigmas.*((abs(z(3)).^2)./Vc);

Uses the vector sigmas, so even if E was a scalar instead of being a vector, the right hand side will be a vector.

Now, if E were a scalar instead of a vector, and if you could make deltac a scalar instead of a vector, then the other variables involved would become scalars and the assignment could work.

ys(2) = (-z(2)./tc)-GLA.*(z(2)-Ntr).*(abs(z(4)).^2)./Vnc;

GLA is a vector constructed in the main routine based on E, so if E were a scalar then GLA would be a scalar.

ys(3) = 1/2*(1-1i)*GN.*(z(1)-Nss).*z(3)+(1/tin.*((z(4)./rR)-z(3)));

Uses the undefined variable Nss. Uses rR which as explained above is a vector because it is built from the vector deltac so if deltac were made into a scalar somehow, then rR would become a scalar and the right hand side would be a scalar (assuming that the undefined Nss is a scalar.)

ys(4) =((-1i.*deltawc-gamat).*z(4))-(P.*gamac.*z(3))+(1/2*(1-1i).*confinec*vg*gN.*(z(2)-Ntr).*z(4));

This used the vectors gamat and gammac from the main routine.

gamai = omega./2/Qi;
gamap = omega./2/Qp;
gamat = omega./2/Qt;

Those use all of the vector omega, so all three are vectors.

gamac = gamat-gamai-gamap;

The right hand side are all vectors, so gamac is a vector.

So ys(4) is a vector.

Beyond defining the undefined n, Nss, and param_rate_eq_fano, there are two approaches you could use:

1. Loop the whole ode45 call over individual values implied by omega, such as by passing in an index and indexing with that were appropriate so that you get 4 x 1 output from the function; Or
2. Reconstruct the ode45 not as a 4 x 1 system but instead as a (4*2001) x 1 system, running all of those 2001 odes at the same time.

[mohammad heydari]

Thank you very much.

actually:

n=3.5;
Nss=(I./(E.*VP))*tc

And I mean param_rate_eq_fano the same constant.I meant calling variables

2.Reconstruct the ode45 not as a 4 x 1 system but instead as a (4*2001) x 1 system, running all of those 2001 odes at the same time.

Do you mean a array should be used?

And what is the number 2001 based on which criterion?

[Walter Roberson]

2001 is length(lambda)

[mohammad heydari]

I tried this too but the problem didn't go away.

thank you.

Answer 2

[Walter Roberson]

The plotting needs to be fixed on this. And it takes a long time!!

 function rate_eq_program
    %%%%%%%%%%%%%%%%%%%%%%%% constants %%%%%%%%%%%%%%%%%%%%%%%%%
    q = 1.6021766208e-19;                       
    c = 2.99792458e8;                                                         
    kk = 1.38e-23;                                                          
    kk_eV = kk/q;                                                            
    T = 300;                                                                  
    h = 6.626068e-34;                                                         
    h_eV = h/q;                                                               
    hbar = h/(2*pi);                        
    hbar_eV = hbar/q;      
    epsilon0 = 8.854187817e-12;  
    P = -1; 
    rL = 1;
    rB = 0.2 ;
    tB = sqrt(1-rB^2);
    Qt =500;
    % Qc =180;
    Qi = 14300;
    Qp = 10000;
    %Qc = 1800;
    L = 5e-6;
    ng = 3.5;
    sigma = 3;
    m = 25;
    A = 0.105e-12; 
    confine = 0.5;
    Vp = 0.576e-18; 
    Vc = confine*A*L;
    gN =5e-20;
    Ntr = 1e24;
    tc = 0.5e-9; 
    %tc = 0.7e-9; 
    alphai = 1000;
    %alpha = 1.5;
    alpha =1;
    confinec = 0.3;
    vg = c/ng;
    vgp = c/6/ng;
    Vpc = 0.22e-18;
    KD = confinec*vgp*alpha*gN/2; 
    KA = -confinec*vgp*gN/2;
    lambda0 = 1554e-9;
    deltalambda = 0.01e-9;
    lambda_start = 1540e-9;
    lambda_end = 1560e-9;
    lambda = lambda_start:deltalambda:lambda_end;
    omega = 2*pi*c./lambda;
    omega0 = omega(lambda==lambda0); %MARK
    E = hbar_eV.*omega;
    gamai = omega./2/Qi;
    gamap = omega./2/Qp;
    gamat = omega./2/Qt;
    %gamac = omega./2/Qc;
    %gamac = omega./2/Qc;
    gamac = gamat-gamai-gamap;
    gama1 = gamac./2;
    gama2 = gama1;
    COS2 = 0.5*gama1./gama1.*tB^2/rB-0.5*tB^2/rB-rB;
    SIN2 = P*tB./(2.*gama1*rB).*sqrt(4.*gama1.*gama1 - tB^2.*(gamac.^2));
    % CS12 = 1/(1i*tB).*(1i*SIN2 + rB);
    GLA = gN*vgp./(E.*Vpc*q);
    GN = confine*vg*gN;
    tin = 2*L/vg;
    M = 0.5*(length(E)-1);
    Tmin = 0;
    Tmax = 5e-9;
    dt = 0.01e-12;
    Time = Tmin:dt:Tmax;
    N0 = zeros(1,4);
    N0(3) = 0.2;
    %I =2e-3;TA I =12e-3;
    I =0.5e-3;
    Vnc= 0.24e-18;
    n=3.5;
    Nss=(I./(E.*Vp))*tc;
    
    NO = length(omega);
    
    %% main code
    %%%%%%%%%%%%%%%%%%%%%%% main code  %%%%%%%%%%%%%%%%%%%%%%%%%%%
    Time = Tmin:dt:Tmax;      
    z0 = repmat([0.2;0.2;0.2;0.2], 1, NO); 
    [t,z] = ode45(@rate_eq_program_1,Time,z0);
    param_rate_eq
    
    %all the plotting is wrong.
    
    figure(1)
    plot(t,z(:,1));   % divided to normalize
    xlabel('time [ns]','FontSize',14);  % size of x label
    ylabel('Arbitrary units','FontSize',14);    % size of y label
    set(gca,'FontSize',14);                 % size of tick marks on both axis
    legend('N', 'Location','SE')  % legend inside the plot
    figure(2)
    plot(t,z(:,2));   % divided to normalize
    xlabel('time [ns]','FontSize',14);  % size of x label
    ylabel('Arbitrary units','FontSize',14);    % size of y label
    set(gca,'FontSize',14);                 % size of tick marks on both axis
    legend('NC', 'Location','SE')  % legend inside the plot
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function %%%%%%%%%%%%%%%%%%%%%%%%
    function ys = rate_eq_program_1(t,z)
    %{
    param_rate_eq_fano
     %}
    z = reshape(z, 4, []);
    
    ys = zeros(size(z));
    
    %rR=(-p*gamma_c)/(1i*(delta_c)+gamma_T-1/2*(1-1i*henry)*conf_NC*V_g*g_n*(z(2)-N_0));
    deltawc=KD.*(z(2,:))-1i*KA.*(z(2,:));
    deltac=omega0+deltawc-omega;
    rR=(-P.*gamac)./(1i*(deltac)+gamat-(1/2*(1-1i)*confinec*vg*gN).*(z(2,:)-Ntr));
    rRS=-P.*gamac./(1i.*(deltac)+gamat);
    %rRS=rB+(-1i*tB-rB).*gamac./(1i.*(deltac)+gamat);
    g=alphai+(1/(2*L)).*log(1./(abs(rRS)).^2);
    A=(2*epsilon0*n*ng)./(hbar.*omega0);
    B=(1+(abs(rRS))).*(1-(abs(rRS)));
    D=((g)-alphai);
    H=c./(omega0.*n);
    F=imag(rRS)./abs(rRS);
    sigmas=A.*((B./D)+(H.*F));
    ys(1,:) = I./(E.*Vp)-(z(1,:)./tc)-GN.*(z(1,:)-Ntr).*sigmas.*((abs(z(3,:)).^2)./Vc);
    ys(2,:) = (-z(2,:)./tc)-GLA.*(z(2,:)-Ntr).*(abs(z(4,:)).^2)./Vnc;
    ys(3,:) = 1/2*(1-1i)*GN.*(z(1,:)-Nss).*z(3,:)+(1/tin.*((z(4,:)./rR)-z(3,:)));
    ys(4,:) =((-1i.*deltawc-gamat).*z(4,:))-(P.*gamac.*z(3,:))+(1/2*(1-1i).*confinec*vg*gN.*(z(2,:)-Ntr).*z(4,:));
    ys = ys(:);
    end
    end

[Walter Roberson]

You have very bad scaling problems. You divide by very small quantities a lot, giving large results, and your GLA is on the order of 1e25 so your ys(2,:) values end up in the 1e24 range. Because the ys(2,:) values are so large, ode45 takes very small steps to try to compute it accurately -- after 9 million function evaluations it was barely up to 1.1e-36 .

You have no hope with this system with ode45. Maybe you would get further with a "stiff" solver such as ode15s .

[mohammad heydari]

Thank you dear Walter.

In fact, given numbers are not approximate numbers, and the numbers are the result of work in laser simulation.I have asked the people who have done this simulation.All of them have done this simulation with ode45.They say that the duration of this simulation is under one minute !! . I really don't know what is going to happen in under a minute !!!

I also got the constant values and initial conditions of this program from someone who simulated the program and reached the final results.In addition, this program you edited is running errors on my system. Does it require a lot of memory?

Best regrads

[Walter Roberson]

I guarantee that none of them used the source code that you gave us.

[mohammad heydari]

I want to put constants together with the formulas. I've worked with a lot of lasers before this program, which took a very long time to simulate. Very powerful systems with a lot of memory needed to be simulated.There are, of course, lasers with very short simulations.This goes back to the values of constant and ode.

[mohammad heydari]

Good time

As I said, I want to put the values here.

In fact, there are other values in the program.It is noteworthy, however, that some indicators are only nominally different from those in the program.Like N0 in program Ntr.

[Walter Roberson]

Your planck constant is not accurate. It is 6.62607015e-34 by definition https://en.wikipedia.org/wiki/Planck_constant but you were using 6.626068e-34 (which I do see in some other places.)

Your Boltzman constant could be more precise. You are using 1.38e-23 when you could be using 1.380649e-23

[Walter Roberson]

I started annotating your constants, but I got bored. You should continue so that people can understand what you are talking about. Make it easy for people to work with your code. Make it easy for you to work with your code, as you will not remember some of these things if you take a long weekend off.

I do not see the Laser cavity cross-section variable ?

%%%%%%%%%%%%%%%%%%%%%%%% constants %%%%%%%%%%%%%%%%%%%%%%%%%
    q = 1.6021766208e-19;       %q = electron charge, coulombs           
    c = 2.99792458e8;           %c = speed of light, m/s                                              
    kk = 1.380649e-23;          %k or $k_B$ Bolzmann constant, J/K              
    kk_eV = kk/q;                                                            
    T = 300;                    %temperature, K                 
    %{
    h = 6.626068e-34;
    %}
    h = 6.62607015e-34;          %planck constant, J*s, definition
    
    h_eV = h/q;                                                               
    hbar = h/(2*pi);                        
    hbar_eV = hbar/q;      
    epsilon0 = 8.854187817e-12;  %\epsilon_0$ = permittivity of vacuum, F/m
    P = -1; 
    rL = 1;            %$R_l$ = Left mirror reflectivity, unit-less
    rB = 0.2 ;
    tB = sqrt(1-rB^2);
    Qt = 500;          %Q_T = Total Q-factor, unit-less
    % Qc =180;
    Qi = 14300;        %Q_i = Intrinsic Q-factor, unit-less
    Qp = 10000;        %Q_p = cross-port Q-factor, unit-less
    %Qc = 1800;
    L = 5e-6;          %L = waveguide cavity length, m
    ng = 3.5;          %$\eta_g$ = group refractive index, unit-less
    sigma = 3;
    m = 25;
    A = 0.105e-12;     %2->0 quenching rate coefficient ??
    confine = 0.5;
    Vp = 0.576e-18; 
    Vc = confine*A*L;
    gN =5e-20;
    Ntr = 1e24;
    tc = 0.5e-9; 
    %tc = 0.7e-9; 
    alphai = 1000;      %$\alpha_i$ = internal loss factor, 1/m
    %alpha = 1.5;
    alpha = 1;          %$\alpha$ = Linewidth enhancement factor, unit-less
    confinec = 0.3;
    vg = c/ng;
    vgp = c/6/ng;
    Vpc = 0.22e-18;
    KD = confinec*vgp*alpha*gN/2; 
    KA = -confinec*vgp*gN/2;
    lambda0 = 1554e-9;        %$\lambda_r$ = Col resonance wavelength
    deltalambda = 0.01e-9;
    lambda_start = 1540e-9;
    lambda_end = 1560e-9;
    lambda = lambda_start:deltalambda:lambda_end;
    omega = 2*pi*c./lambda;
    omega0 = omega(lambda==lambda0); %MARK

    omega = omega0;
    
    E = hbar_eV.*omega;
    gamai = omega./2/Qi;
    gamap = omega./2/Qp;
    gamat = omega./2/Qt;
    %gamac = omega./2/Qc;
    %gamac = omega./2/Qc;
    gamac = gamat-gamai-gamap;
    gama1 = gamac./2;
    gama2 = gama1;
    COS2 = 0.5*gama1./gama1.*tB^2/rB-0.5*tB^2/rB-rB;
    SIN2 = P*tB./(2.*gama1*rB).*sqrt(4.*gama1.*gama1 - tB^2.*(gamac.^2));
    % CS12 = 1/(1i*tB).*(1i*SIN2 + rB);
    GLA = gN*vgp./(E.*Vpc*q);
    GN = confine*vg*gN;
    tin = 2*L/vg;
    M = 0.5*(length(E)-1);
    Tmin = 0;
    Tmax = 5e-9;
    dt = 0.01e-12;
    Time = Tmin:dt:Tmax;
    N0 = zeros(1,4);
    N0(3) = 0.2;
    %I =2e-3;TA I =12e-3;
    I =0.5e-3;
    Vnc= 0.24e-18;
    n=3.5;
    Nss=(I./(E.*Vp))*tc;

[Walter Roberson]

I ran a modified version of the code that restricts omega to one value so that the variable sizes all match up in the function. Almost immediately it started needing to use a timestep of approximately 1e-37 to try to deal with how steep the variables are -- since you do, after-all, go from a boundary condition of 0.2 to an output value above 1e24. With a time-step of about 1e-37, you would take a very long time to finish.

When ode45 needs very fine time steps, it is often an indication that the system is "stiff" and needs to use a stiff solver, ode15s or ode23s.

When I use ode23s, it gives a large number of warnings about matrix being singular, but makes it through the complete time span. The z outputs are all complex-valued, and it is difficult to see how any of them could be useful.

When I use ode15s, it does not give any warnings about matrix being singular, but it stops at about time 1e-10 saying that it has probably encountered a singularity and cannot meet integration tolerances. The partial output it produces when it stops has z outputs that are all complex-valued and it is difficult to see how any of them could be useful.

[mohammad heydari]

I thank you and I apologize for the delay in responding and not giving you enough explanation for the same reasons you said.I understand your point, but there are some ambiguities in this regard.As far as you can see from the reviews you have done, the main problems are with the parameters and this can affect the output performance.Of course, I believe that the parameters are reasonably accurate.Because I've looked at these numbers in several authoritative articles and can almost assure you that the parameters have been entered correctly, I've shown you some of them.Can't we get the right output with the changes in the matrix Z and ys?

I say this because ode15s does not give us an accurate approximation of how the laser works.I would like to focus more on the main structure of the program, which includes commands ode45.I'm pretty sure that by changing the program structure, you can get the right outputs.

In addition

A= Laser cavity cross-section

Best regards

[Walter Roberson]

Sorry, I do not appear to be able to communicate successfully with you on this matter. Perhaps some other volunteer will have more success. However, in my experience, most of the other volunteers are much stricter about requiring commented code than I am.

[mohammad heydari]

thank you very much for your help.

Of course, the disagreement is normal in these cases. And I'll try to include more explanations along with the formulas over the coming hours to give a better description of the situation.

In any case, I sincerely thank you