Solving for a system

This code works well when ax takes only a single value. I wanted to see the dependence of Lc as the value of ax changes. But when ax takes a range of values as given here I get an error. What changes do I have to make here?

lambda=1.55e-6;

%Field in vertical direction
ax=linspace(1e-6,2e-6,10);
nw=1.5;
ns=1.498;
% Find value of V0
V0=2*pi/lambda*ax*(nw^2-ns^2)^.5;
syms b;
ri=(nw/ns)^2;
eqn0= V0*(1-b)^.5==2*atan(ri*(b/(1-b))^.5);
b0=double(vpasolve(eqn0,b));
neff0=(b0*(nw^2-ns^2)+ns^2)^.5;
kappax=2*pi/lambda*(nw^2-neff0^2)^.5;
gammasx=2*pi/lambda*(neff0^2-ns^2)^.5;
CMeqn_lx=1/gammasx+(ax+sin(kappax*ax)/kappax)/2/(cos(kappax*ax/2)^2);
%Field in horizontal
ay=8e-6;
gapy=2e-6;
V=2*pi/lambda*ay*(neff0^2-ns^2)^.5;
ri=1;
eqn= V*(1-b)^.5==2*atan(ri*(b/(1-b))^.5);
b1=double(vpasolve(eqn,b));
neff=(b1*(neff0^2-ns^2)+ns^2)^.5;
beta=2*pi/lambda*neff;
kappay=2*pi/lambda*(neff0^2-neff^2)^.5;
gammasy=2*pi/lambda*(neff^2-ns^2)^.5;
%Coupling Mode Equation
CMeqn_ly=1/gammasy+(ay+sin(kappay*ay)/kappay)/2/(cos(kappay*ay/2)^2);
CMeqn_r=4*pi/lambda/beta*3e8*(4*pi*1e-7);
C2=CMeqn_r/CMeqn_lx/CMeqn_ly;
C=C2^.5
intgx=(ax+sin(kappax*ax)/kappax)/2/(cos(kappax*ax/2)^2);
intgy=(exp(-gammasy*gapy)*cos(kappay*ay/2)-exp(-gammasy*(ay+gapy))*cos(kappay*ay/2))/gammasy+(exp(-gammasy*gapy)*sin(kappay*ay/2)+exp(-gammasy*(ay+gapy))*sin(kappay*ay/2))*kappay/(gammasy^2);
intgy=intgy/(1+(kappay/gammasy)^2)/cos(kappay*ay/2);
Coupling=2*pi*3e8/lambda/4*(8.85418782e-12)*(nw^2-ns^2)*C2*intgx*intgy
Lc=pi/2/Coupling;

Answer 1

[Star Strider]

The easiest way is to loop through the values of ‘ax’:

lambda=1.55e-6;
%Field in vertical direction
axv=linspace(1e-6,2e-6,10);
nw=1.5;
ns=1.498;
for k = 1:numel(axv)
    ax = axv(k);
    % Find value of V0
    V0=2*pi/lambda*ax*(nw^2-ns^2)^.5;
    syms b;
    ri=(nw/ns)^2;
    eqn0= V0*(1-b)^.5==2*atan(ri*(b/(1-b))^.5);
    b0=double(vpasolve(eqn0,b));
    neff0=(b0*(nw^2-ns^2)+ns^2)^.5;
    kappax=2*pi/lambda*(nw^2-neff0^2)^.5;
    gammasx=2*pi/lambda*(neff0^2-ns^2)^.5;
    CMeqn_lx=1/gammasx+(ax+sin(kappax*ax)/kappax)/2/(cos(kappax*ax/2)^2);
    %Field in horizontal
    ay=8e-6;
    gapy=2e-6;
    V=2*pi/lambda*ay*(neff0^2-ns^2)^.5;
    ri=1;
    eqn= V*(1-b)^.5==2*atan(ri*(b/(1-b))^.5);
    b1=double(vpasolve(eqn,b));
    neff=(b1*(neff0^2-ns^2)+ns^2)^.5;
    beta=2*pi/lambda*neff;
    kappay=2*pi/lambda*(neff0^2-neff^2)^.5;
    gammasy=2*pi/lambda*(neff^2-ns^2)^.5;
    %Coupling Mode Equation
    CMeqn_ly=1/gammasy+(ay+sin(kappay*ay)/kappay)/2/(cos(kappay*ay/2)^2);
    CMeqn_r=4*pi/lambda/beta*3e8*(4*pi*1e-7);
    C2=CMeqn_r/CMeqn_lx/CMeqn_ly;
    C=C2^.5;
    intgx=(ax+sin(kappax*ax)/kappax)/2/(cos(kappax*ax/2)^2);
    intgy=(exp(-gammasy*gapy)*cos(kappay*ay/2)-exp(-gammasy*(ay+gapy))*cos(kappay*ay/2))/gammasy+(exp(-gammasy*gapy)*sin(kappay*ay/2)+exp(-gammasy*(ay+gapy))*sin(kappay*ay/2))*kappay/(gammasy^2);
    intgy=intgy/(1+(kappay/gammasy)^2)/cos(kappay*ay/2);
    Coupling=2*pi*3e8/lambda/4*(8.85418782e-12)*(nw^2-ns^2)*C2*intgx*intgy;
    Lc(k)=pi/2/Coupling;
end

figure
plot(axv, Lc)
grid
xlabel('ax')
ylabel('Lc')