function neff=step_index
a=10e-6;
n1=1.44;
n2=1.4;
lam=1.55e-6;

k0=2*pi/lam;
beta=linspace(n2*k0,n1*k0,1000);

T=TMeq(beta,a,n1,n2,k0);

intervals=(T<.1)&(T>-.1);
istart=find(diff(intervals)>0);
istop=find(diff(intervals)<0);
X0=[beta(istart); beta(istop)]';

[nzeros,scrap]=size(X0);

for i=1:nzeros
    neff(i)=fzero(@(x) TMeq(x,a,n1,n2,k0),X0(i,:))/k0;
end

neff=neff(end:-1:1);

function y=TMeq(beta,a,n1,n2,k0)

q=sqrt(beta.^2-n2^2*k0^2);
h=sqrt(n1^2*k0^2-beta.^2);

T1=real(besselj(1,h*a)./(h*a.*besselj(0,h*a)));
T2=real(n2^2*besselk(1,q*a)./(q*a*n1^2.*besselk(0,q*a)));

y=T1+T2;