rm = 6e-32
rl = 20e-9
vo = 0.8
h = 6.626e-34
pi=3.141592654
hb = h/2./pi
q = 1.602e-19

i = 1

for E = 0.001:0.001:(vo-0.001)
    rhs(i) = sqrt(vo/E-1)
    rlhs_factor = rl*sqrt(2.*rm*q/hb*E/hb)/2.
    rlhs1(i) = tan(rlhs_factor)
    rlhs2(i) = -cot(rlhs_factor)    
    xvec(i)=E
    i = i + 1
end

plot(xvec,rhs,xvec,rlhs2);

for  i=1:1:8
    en =hb/q*hb/2./rm*(real(i)*pi/rl)^2
    i,en
end