clear all, close all, clc
n=1000
rmax=10; r0=-rmax;h=2*rmax/(n); 
emin=1.1; emax=1.5; he=.1;p=1; 
for e=emin:he:emax
x=[1;0];
for j=n:-1:0           
rd=r0+j*h;k1=sqrt(2*(e-sig(rd))); k2=sqrt(2*(e-sig(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)=rexas(e);
ener(p)=e;p=p+1;
end
[ener',rcal',rext']
%plot(ener,rcal,'-')