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

error

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.

Best Answer

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.
Related Question