clear all, close all, clc
n=100
rmax=10; r0=-rmax;h=2*rmax/(n); 
emin=.1; emax=3; he=.011;p=1; 
for e=emin:he:emax
x=[1;0];
for j=n:-1:0           
rd=r0+j*h;k1=sqrt(2*(e-bar(rd))); k2=sqrt(2*(e-bar(rd+h)));
P2=(-1/(2*k1))*[-k1*exp(-i*k1*rd), -exp(-i*k1*rd);...
-k1*exp(i*k1*rd), exp(i*k1*rd)];
P1=[exp(i*k2*rd), exp(-i*k2*rd);...
k2*exp(i*k2*rd), -k2*exp(-i*k2*rd)];
P=P2*P1;x=P*x;
end
rcal(p)=(x(2)*x(2)')/(x(1)*x(1)');
rext(p)=rexa(e);
ener(p)=e;p=p+1;
end
[ener',rcal',rext'];
plot(ener,1-rcal,'-',ener,1-rext,':')