# 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_eqfigure(1)plot(t,z(:,1));   % divided to normalizexlabel('time [ns]','FontSize',14);  % size of x labelylabel('Arbitrary units','FontSize',14);    % size of y labelset(gca,'FontSize',14);                 % size of tick marks on both axislegend('N', 'Location','SE')  % legend inside the plotfigure(2)plot(t,z(:,2));   % divided to normalizexlabel('time [ns]','FontSize',14);  % size of x labelylabel('Arbitrary units','FontSize',14);    % size of y labelset(gca,'FontSize',14);                 % size of tick marks on both axislegend('NC', 'Location','SE')  % legend inside the plot%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function %%%%%%%%%%%%%%%%%%%%%%%%function ys = rate_eq_program_1(t,z)param_rate_eq_fanoys=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.

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.