n=linspace(1.43,1.50,10);
function y=f(x,na,da)
ea=na.^2;
n will eventually be 1.50 because of the linspace(), so ea = na.^2 will eventually be 2.25
ka=k0*sqrt(ea-x^2);
Suppose that x happens to be 1.5... perhaps because
p1 = 1.5;
then ea-x^2 would be 1.50^2 - 1.5^2 --> 0 so ka would become 0
rac=(ka-kc)/(ka+kc);
ka is 0, so rac is -kc/kc which will be -1 unless kc is exactly 0.
f=(1-rac*exp(2*1i*da*ka))/(1+rac*exp(2*1i*da*ka));
Remember that ka is 0, so 2*1i*da*ka is 0 unless da is infinite or nan. exp(0) is 1. Remember that rac is -1 because it is -kc/kc when ka is 0. So you now have -1*1 which is 1. And then you have 1+(-1) which is 0. So you have a division by 0. Which gives you NaN, which poisons the rest of your calculation.
P(j)=real(p);
You need to double-check that. You created P as a 2D array
