from running code find diff plot vs variable

function mm1
clear all
close all

n=linspace(1.43,1.50,10);
D1= linspace(1e-9,3e-7,15);
P=zeros(numel(n),numel(D1));
for j=1:numel(D1)
    for k=1: numel(n)  
        da = D1(j);
        na=n(k);
        p0 = 1;
        p1 = 1.5;
        p2 = 2;
        TOL = 10^-8;
        N0 = 100; format long
        h1 = p1 - p0;
        h2 = p2 - p1;
        DELTA1 = (f(p1,na,da) - f(p0,na,da ))/h1;
        DELTA2 = (f(p2,na,da ) - f(p1,na,da ))/h2;
        d = (DELTA2 - DELTA1)/(h2 + h1);
        i=3;
        while i <= N0 
            
            b = DELTA2 + h2*d;
            D = (b^2 - 4*f(p2,na,da )*d)^(1/2);
            
            if abs(b-D) < abs(b+D)
                E = b + D;
            else
                E = b - D;
            end
            
            h = -2*f(p2,na,da)/E;
            p = p2 + h;
            
            if abs(h) < TOL      
                disp(p) 
                break
            end
            
            p0 = p1;
            p1 = p2;
            p2 = p;
            h1 = p1 - p0;
            h2 = p2 - p1;
            DELTA1 = (f(p1,na,da ) - f(p0,na,da ))/h1;
            DELTA2 = (f(p2,na,da ) - f(p1,na,da ))/h2;
            d = (DELTA2 - DELTA1)/(h2 + h1);
            i=i+1;
            
        end
    end
    %pv(j)=p;
    P(j)=real(p);
    %R(j)=real(r);
end
plot(D1,P)
end

function  y=f(x,na,da)

l=6330*(10)^(-10);

k0=2*pi/l;
nc=1.33;
ec=nc^2;
%na=1.43;
ea=na.^2;
nf=1.59;
ef=nf^2;
nm=0.064+1i*4;
em=nm^2;
ns=1.495;
es=ns^2;
df=0.5e-6;
dm=0.03e-6;
%kx=(2*pi/l)*np*sin(a);
kc=k0*sqrt(x^2-ec);
ka=k0*sqrt(ea-x^2);
kf=k0*sqrt(ef-x^2);
km=k0*sqrt(em-x^2);
ks=k0*sqrt(x^2-es);
rfm=(kf-km)/(kf+km);
rms=(km-ks)/(ks+km);
rfa=(kf-ka)/(kf+ka);
rac=(ka-kc)/(ka+kc);
a=(1-rfm)/(1+rfm);
b=(1-rms*exp(2*1i*dm*km))/(1+rms*exp(2*1i*dm*km));
c=(1-rfa)/(1+rfa);
f=(1-rac*exp(2*1i*da*ka))/(1+rac*exp(2*1i*da*ka));
pfms=2*atan(1i*a*b);
pfac=2*atan(1i*c*f);
y=2*kf*df-pfms-pfac;






end

pl help to plot between dy/dna vs da

Answer 1

[Walter Roberson]

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