Dear all,
Please help me with the next problem.I'm beginner!
I made the following algorithm:

c=3e8;//m/s
h=6.6e-34;//J*s
lambda=200e-9;//m
E=h*c/lambda;//J
k0=(2*E)/(h*c);//m^-1
l=10;
n=1.33;
rd=20*10^-6;
r=1*10^-9:200:10*rd;//m
for i=1:length(r)
if r(i)<=rd then
Ep(i)=(E*(n^2-1))+(E*l*(l+1)/(k0^2*r(i).^2));//J
elseif r>rd then
Ep(i)=(E*(n^2-1))+(E*l*(l+1)/(k0^2*r(i).^2));
end
end
clf
plot2d(r,Ep,style=1);
xlabel("$rd\le r\le 10*rd$","fontsize",4,"color","red");
ylabel("Ep","fontsize",4,"color","red");
title("Ep=f(r)","color","red","fontsize",4);
p = get("hdl");
p.children.mark_mode = "on";
p.children.mark_style = 4;
p.children.thickness = 1;
p.children.mark_foreground = 5;

I must obtain a curve and I don't know where is the error.For the moment my program calculates only a value.Please tell me if you see the error.
Thanks!