clear all, clc, close all
n=1000
rmax=10; r0=-rmax;h=2*rmax/(n); 
emin=.1; emax=4.5; he=.0123; 
a1=1:0.5:2.5;
cc=0;
for a=a1
p=1;
for e=emin:he:emax
x=[1;0];
rr=0;
for j=n:-1:0           
rd=r0+j*h;k1=sqrt(2*(e-eckV3(a,rd))); k2=sqrt(2*(e-eckV3(a,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;

rr=rr+1;
R(rr)=rd;
end
rcal(p)=(x(2)*x(2)')/(x(1)*x(1)');
rext(p)=rexaeck(a,e);
ener(p)=e;p=p+1;
end
cc=cc+1;
[ener',rcal']
figure(1)
plot(ener,rcal,'-',ener(1:6:p),rext(1:6:p),'o'),hold on,
txt = ['a = ',num2str(a)];
ind=[92 143 194 245];
text(ener(ind(cc)), rcal(ind(cc)), txt,'horizontalalignment','left')
xlabel('E')
ylabel('R')
%print -deps -r900 RvsEeck

figure(2)
plot(R,eckV3(a,R),'-'),hold on,
txt = ['a = ',num2str(a)];
text(R(110), eckV3(a,R(110)), txt,'verticalalignment','bottom')
xlabel('x')
ylabel('E_p(x)')
%print -deps -r900 EckPot
end