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
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 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
g is a vector so D will be a vector.
Use the undefined variable n. If we assume that n is a scalar, then H will be a scalar.
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:
- 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
- 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.
