
n1=1;
n2=1.5;
function n=indice2(z)
    n=(n1+z*(n2-n1))^2;
endfunction
lambda=10;
            

k=2*%pi/lambda;
function [deriv]=systeme(z,y)
    deriv(1)=k*y(4),
    deriv(2)=-k*y(3),
    deriv(3)=k*indice2(z)*y(2),
    deriv(4)=-k*indice2(z)*y(1),
endfunction
            

z=[0:0.01:1];
tolA=1d-6;
tolR=1d-10;
V1=ode([1;0;0;0],0,z,tolR,tolA,systeme);
V2=ode([0;0;1;0],0,z,tolR,tolA,systeme);
p=length(z);
M=[V2(3,p)+%i*V2(4,p), -V2(1,p)-%i*V2(2,p); -V1(3,p)-%i*V1(4,p), V1(1,p)+%i*V1(2,p)]
            

function r=reflexion(M,ni,nt)
    r=((M(1,1)+M(1,2)*nt)*ni-(M(2,1)+M(2,2)*nt))/((M(1,1)+M(1,2)*nt)*ni+(M(2,1)+M(2,2)*nt))
endfunction
function t=transmission(M,ni,nt)
    t=2*ni/((M(1,1)+M(1,2)*nt)*ni+(M(2,1)+M(2,2)*nt))
endfunction
            
r=abs(reflexion(M,n1,n2))^2
