clear all
h = 6.62607e-34/(2*pi); %Js
m = 9.109383e-31; %kg

a = 1e-9; %m
V0 = 5*1.602e-19; %J

E = linspace(-V0,0,10001);

phi = sqrt(2*m*a^2*(E+V0)/h^2);
display([ 'phi_{max} = ' num2str(max(phi)/(pi),'%2.2f') ' pi' ]) %Max fas

f1 = sqrt(-E./(E+V0));
f2 = tan(phi);

figure(1)
plot(E/1.602e-19,[f1; f2],'.')
axis([min(E)/1.602e-19 0 -10 10])
ylabel('tan(\phi) and (-E/(E+V_0))^{1/2}')
xlabel('E (eV)')
