k0=(2*pi/0.6328)*1e6;
t2=1.5e-6;
n1=1.512;n2=1.521;n3=4.1-1i*0.211;
n4=1;
m=0;
x=2;
t3=linspace(1e-9,1e-6,20)
k1=k0*sqrt(n1^2-x.^2);k2=k0*sqrt(n2^2-x.^2);
k3=k0*sqrt(n3^2-x^2);
k4=k0*sqrt(n4^2-x.^2);
tol = 1e-12;
n = 1;
while (n<=200)
x_new =-(k2).*t2+atan(k1./1i*k2)+atan((k3./k2).*tan(atan(k4./1i*k2)-k3.*t3))+m*pi;
if abs(x_new-x) < tol, break; end
x =x_new;
k1=k0*sqrt(n1^2-x.^2);
k2=k0*sqrt(n2^2-x.^2);
k3=k0*sqrt(n3^2-x.^2);
k4=k0*sqrt(n4^2-x.^2);
n=n+1;
end
disp(x_new)
plot(t3,real(x_new))
edit program