# set the FDE solver as active
switchtolayout;
setactivesolver("FDE");
# find mode of the single input waveguide
setnamed("FDE","x",-3e-6);
findmodes;
cleardcard("E0");
copydcard("mode1","E0");
# propagate in MMI along x direction and collect E intensity
L = linspace(0,6e-6,200); # calculate the fields at 200 positions over 6 um distance
y = getdata("mode1","y");
z = getdata("mode1","z");
E2 = matrix(length(L),length(y)); # initialize matrix to store E intensity
# move solver to cross section of MMI and find modes
switchtolayout;
setnamed("FDE","x",0);
setanalysis("number of trial modes",20);
findmodes;
# minimum and maximum effective index of modes to include in overlap calculations when propagating
nmin = 0.1;
nmax = 4;
for(i=1:length(L)) { # loop over positions along x
outmode = propagate("E0",L(i),nmin,nmax);
E2_temp = getelectric(outmode);
cleardcard(outmode);
E2(i,1:length(y)) = E2_temp(1,1:length(y),find(z,0.1e-6));
}
# plot image of propagated E field intensity
image(L*1e6,y*1e6,E2,"L (microns)","y (microns)","|E|^2 vs L");