Determining Coupling coefficent value for mode coupling in a dual core fiber solving characterstics equation

The graph obtained is not as expected. Can someone have a look on my code and figure out mistakes in it?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lambda=(1:0.01:2); %unit: um

%constants and variables

pi=3.14159; %Greek letter pi

c=2.99792458e14; % speed of light (unit: um/s)

a=1.95; % radius of core (unit: um)

d=2.52*4; % distance of separation between cores (unit: um)

d1=1.1; % air hole diameter (unit:um)

P=2.52; %pitch or hole-hole spacing

ko=2*pi./lambda;

dbar=d/a;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Calculate refractive index of core material (n1) using Sellmeier's equation

G1=0.6961663; % Sellmeier's coefficients for pure silica

G2=0.4079426;

G3=0.8974794;

lambda1=.0684043; %wavelength for pure silica

lambda2=0.1162414;

lambda3=9.896161;

npower2oflambda=1+(G1*lambda.^2./(lambda.^2-power(lambda1,2)))+(G2*lambda.^2./(lambda.^2-power(lambda2,2)))+...

(G3*lambda.^2./(lambda.^2-power(lambda3,2)));

noflambda=sqrt(npower2oflambda);

n1=noflambda;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%calculate n2 in terms of n1

NA=0.20;

n2=sqrt(n1.^2-NA^2);

% calculate V-number

V=ko.*a.*(n1.^2-n2.^2 ).^(1/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%SOLVING BESSEL FUNCTION EIGEN VALUE EQUATION

syms U W

k1=(ko.*n1); %maxm value of beta

k2=(ko.*n2); %minimum value of beta

% beta=linspace(k2,k1,101);

beta=real(k2):real(k1-k2)/1e2:real(k1);

% n_eff=real(n2):real(n1-n2)/1e2:real(n1);

% beta=ko.*n_eff;

U=a.*sqrt(ko.^2.*n1.^2-beta.^2); %transverse mode paramter (phase constant)

W=a.*sqrt(beta.^2-ko.^2.*n2.^2); %transverse mode parameter (attenuation coefficient)

% h=U.*besselj(1,U)./besselj(0,U);

% m=W.*besselk(1,W)./besselk(0,W);

%solution=solve(U.*besselj(1,U)./besselj(0,U)==W.*besselk(1,W)./besselk(0,W), W);

%solution=double(solution)

figure(1)

plot(beta,h,'g','linewidth',2)

hold on

plot(beta,m,'r','linewidth',2)

legend('J_1(U)/U*J_0(U)','K_1(W)/W*K_0(W)')

xlabel('\delta\beta(um)')

grid on

hold off

betaroot=9.034; % first root from the plot 1

U1=a*sqrt(ko.^2.*n1.^2-betaroot.^2); %transverse mode parameter (phase constant)

W1=a*sqrt(betaroot.^2-ko.^2.*n2.^2);

%Uroot=imag(U1);

%Wroot=abs(W1);

K_0=besselk(0,W1.*dbar); %Zeroth order modified bessel function of 2nd kind

K_1=besselk(1,W1); %first order modified bessel function of 2nd kind

% %calculate coupling coefficient

R=sqrt(2*del)/a.*((U1.*U1)./(V.*V.*V)); % 1st part of coupling coefficient equation

Q=K_01./(K_1.*K_1); % 2nd part of coupling coefficient equation

q=abs(Q);

CC=R.*q;%Q;

figure(2)

plot(lambda,CC)

xlabel('wavelength (um)'); ylabel('coupling coeff. (/m)')

title('Coupling coeff. Solving characterstics equation');