% Constants specification
m = 6e-32
q = 1.602e-19
h = 6.625e-34
pi = 3.141592654
hb = h/pi/2
ri = 1i

% Potential well parameters specification
V1 = 0
V2 = 0.5
V3 = 0
V4 = 0.5
V5 = 0
L1 = 4e-9
L2 = 4e-9
W =5e-9

i = 1
for E = 0.01:0.001:1

%   Calculate k-vectors
    k1 = sqrt(2.*m/hb*q/hb*(E-V1))
    k2 = sqrt(2.*m/hb*q/hb*(E-V2))
    k3 = sqrt(2.*m/hb*q/hb*(E-V3))
    k4 = sqrt(2.*m/hb*q/hb*(E-V4))
    k5 = sqrt(2.*m/hb*q/hb*(E-V5))

%   Calculate Matrices
    M1 = ([(k1+k2)/2/k1 (k1-k2)/2/k1;(k1-k2)/2/k1 (k1+k2)/2/k1]);
    factor = exp(-ri*k2*L1)
    M2 = ([factor 0; 0 1/factor]);
    M3 = ([(k2+k3)/2/k2 (k2-k3)/2/k2;(k2-k3)/2/k2 (k2+k3)/2/k2]);
    factor = exp(-ri*k3*W)
    M4 = ([factor 0; 0 1/factor]);
    M5 = ([(k3+k4)/2/k3 (k3-k4)/2/k3;(k3-k4)/2/k3 (k3+k4)/2/k3]);
    factor = exp(-ri*k4*L2)
    M6 = ([factor 0; 0 1/factor]);
    M7 = ([(k4+k5)/2/k4 (k4-k5)/2/k4;(k4-k5)/2/k4 (k4+k5)/2/k4]);
    MT = M1*M2*M3*M4*M5*M6*M7
    trans(i) = 1/abs(MT(1,1))^2
    refle(i) = abs(MT(2,1)/MT(1,1))^2
    energy(i) = E
    i = i+1
end

figure(1);
semilogy(energy,trans);
figure(2)
semilogy(energy,refle);

clear all;