clear all
m_0 = 9.10E-31; % in Kg
m_eff = 0.07*m_0;

m1_eff = m_eff;
m2_eff = m_eff;
h = 4.135E-15; % planks in eV
%h = 6.626E-31;
h_bar = h/(2*pi);
h_bar_sq = h_bar^2;
e = 1.6E-19;
V = 1;
a = h_bar*10/sqrt(2*m_eff*V);
%a = 100E-9;
E_range = 0:0.01:4;
count = 1;
for E = E_range
k0 = sqrt(2*m_eff*E)/h_bar_sq;
if E < V
    k1 = sqrt(2*m_eff*(V-E))/h_bar_sq;
elseif E > V
    k1 = sqrt(2*m_eff*(E-V))/h_bar_sq;
else
    k1 = k0;
end
Y = [1;0;k0;0];
M(:,1) = [-1;0;k0;0];
M(:,2) = [1;-exp(1i*k1*a);k1;-k1*exp(1i*k1*a)];
M(:,3) = [1;-exp(-1i*k1*a);-k1;k1*exp(-1i*k1*a)];
M(:,4) = [0;exp(1i*k0*a);0;k0*exp(1i*k0*a)];
%M = transpose(M);
X(:,count) = M\Y;
if E<V
    T(count) = 1/(1 + ((V^2*(sinh(k1*a))^2)/(4*E*(V-E))));
elseif E >V
    T(count) = 1/(1 + ((V^2*(sin(k1*a))^2)/(4*E*(E-V))));
else
    T(count) = 1/(1 + (m_eff*a^2*V/2*h_bar_sq));
end
    t(count) = abs(4*k0*k1*exp(-1i*a*(k0 - k1))/((k0 + k1)^2 - exp(2i*a*k1)*(k0 - k1)^2))^2;
    count = count +1;
end


plot(E_range,(abs(X(4,:))))
hold on
plot(E_range,t)