function Guided_Modes(d,lamda)

n1 = 2.3278;
n2 = 2.2028;
n3 = 1.00;
k0 = 2*pi/lamda;
k1 = k0*n1;
aE = (n2^2-n3^2)/(n1^2-n2^2);
aM = (n1^4/n3^4)*(n2^2-n3^2)/(n1^2-n2^2);
V = k0*d*sqrt(n1^2-n2^2);

hdTE = [];
hdTM = [];
BgTE = [];
BgTM = [];
x = 0;
    while x+3 < V;
        syms h1d
        a2d = sqrt(V^2-h1d^2);
        a3d = sqrt((V^2)*(1+aE)-h1d^2);
        a3dm = sqrt((V^2)*(1+aM)-h1d^2);
        f = tan(h1d);
        g = (h1d*(a2d+a3d))/(h1d^2-a2d*a3d);
        q = (h1d/n1^2*(a2d/n2^2+a3dm/n3^2))/(h1d^2/n1^4-(a2d*a3dm)/(n2^2*n3^2));
        
        TE = f-g;
        TM = f-q;
        TEmode = vpasolve(TE == 0,h1d, [0.1+x,3+x]);
        hdTE = [hdTE;TEmode];
        TMmode = vpasolve(TM == 0,h1d, [0.1+x,3+x]);
        hdTM = [hdTM;TMmode];

        a2d = sqrt(V^2-TEmode^2);
        a2dm = sqrt(V^2-TMmode^2);
        a3d = sqrt((V^2)*(1+aE)-TEmode^2);
        a3dm = sqrt((V^2)*(1+aM)-TMmode^2);
        
        BTE = sqrt(k1^2-TEmode^2/d^2);
        BgTE = [BgTE;BTE];
        phiTE = atan((TEmode*(a2d-a3d)/(TEmode^2+a2d*a3d)))/2
        neffTE = n1*BgTE/k1;

        BTM = sqrt(k1^2-TMmode^2/d^2);
        BgTM = [BgTM;BTM];
        phiTM = atan((TMmode/n1^2*(a2dm/n2^2-a3dm/n3^2))/(TMmode^2/n1^4+(a2dm*a3dm)/(n2^2*n3^2)))/2;
        neffTM = n1*BgTM/k1;
        x = x + 3;
    end;
 hdTE
 BgTE
 hdTM;
 BgTM;
end