numerical Instabilities for bessel functions

I write the following code but i used the following parametrs and i did not get any numerical instabilities:

RI = 34.5; % Inner radius (m)

ratio = 47.5/34.5; % Ratio R_E / R_I

RE = ratio*RI; % outer radius (m)

h = 200/34.5 * RI; % Water depth (m)

d = 38/34.5 * RI; % Interface depth (m)

b = h-d;

But when i switch the parameters with this then i get numerical instabilities

RI = 0.4; % Inner radius (m)

RE = 0.8; % Outer radius (m)

d = 0.3; % Draft (m)

h = 1.0; % Water depth (m)

b = h-d; %from cylinder bottom to seabed (m)

the following is code

clear all;
close all;
clc;
tic;
% diary('iter_log.txt'); 
% diary on
% fprintf('Run started: %s\n', datestr(now));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%          Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 10;                    % Number of N modes
M = 10;                    % Number of M modes
Q = 10;                    % Number of Q modes
RI = 0.4;        % Inner radius (m)
RE = 0.8;        % Outer radius (m)
d = 0.3;         % Draft (m)
h = 1.0;         % Water depth (m)              % Interface depth (m)
b = h-d;                   % Distance from cylinder bottom to seabed (m)
g = 9.81;       % Gravity (m/s^2)
tau =0.2; % porosity ratio
gamma = 1.0;               % discharge coefficient
               
b1    = (1 - tau) / (2 * gamma * tau^2);   % nonlinear  coeff

tol      = 1e-4;                              % convergence tolerance
max_iter = 100;                                % max  iterations
l = 0;                     % Azimuthal order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === Compute Bessel roots and scaled chi values ===
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
roots = bessel0j(l, Q);    % Zeros of J_l(x)
chi = roots ./ RI;         % chi_q^l = roots scaled by radius
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    =====     Define Range for k_0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % X_kRE = linspace(0.0025, 0.16,100);         %%%%%% this is now k_0 
X_kRE = linspace(0.05, 4.5, 40);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE));  % Preallocate
omega = zeros(size(X_kRE));  %%%% omega preallocate
% iters_used = zeros(length(X_kRE),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for idx = 1:length(X_kRE)
     wk = X_kRE(idx);                                % dimensional wavenumber  
        omega(idx) = sqrt(g * wk * tanh(wk * h));   % dispersion relation
        % fprintf('\n--- idx %3d (T=%6.3f s) ---\n', idx, 2*pi/omega(idx));
% drawnow;
%%%%%%%%%%%%% Compute group velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Cg(idx) = (g * tanh(wk * h) + g * wk * h * (sech(wk * h))^2) ...
              * omega(idx) / (2 * g * wk * tanh(wk * h));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% =======Compute a1 based on tau 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    a1(idx) = 0.93 * (1 - tau) * Cg(idx);
    % dissip(idx) = a1(idx)/omega(idx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %=============== Find derivative of Z0 at z = -d 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  Z0d = (cosh(wk * h) * wk * sinh(wk * b)) / (2 * wk * h + sinh(2 * wk * h));
  fun_alpha = @(x) omega(idx)^2 + g*x.*tan(x*h);     %dispersion equation for evanescent modes
  opts_lsq = optimoptions('lsqnonlin','Display','off');  % ← add this
 for n = 1:N
    if n == 1
        k(n) = -1i * wk;
    else

        x0_left  = (2*n - 3) * pi / (2*h);
        x0_right = (2*n - 1) * pi / (2*h);
        guess = (x0_left + x0_right)/2;

        % Use lsqnonlin for better convergence
         k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right,opts_lsq);
         %%%%%%%%%%%%% Derivative of Z_n at z = -d %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Znd(n) = -k(n) * (cos(k(n)*h) * sin(k(n)*b)) / (2*k(n)*h + sin(2*k(n)*h));
       
    end
end
%% finding the root \lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;

 end
% %%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  A = zeros(M,N);
   A(1,1) = (cosh(wk*h)*sinh(wk*b))/(wk*b*(2*wk*h+sinh(2*wk*h)))*(besselh(l, wk*RE)) /(dbesselh(l, wk*RE)); % Set A_00 
 for j = 2:N
    A(1,j) =(cos(k(j)*h)*sin(k(j)*b) )/ (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))*( besselk(l, k(j)*RE) )...
       /(dbesselk(l, k(j)*RE)); % Set A_0n

 end


for i = 2:M
  A(i,1) =(2 * cosh(wk*h) * (-1)^(i-1) * wk*sinh(wk*b)) / (b * (2*wk*h + sinh(2*wk*h))*(wk^2 + m(i)^2))...
      * (besselh(l, wk*RE))/ (dbesselh(l, wk*RE)); % Set A_m0 
end
for i = 2:M
    for j = 2:N
        A(i,j) = (2*cos(k(j)*h)*(-1)^(i-1) *k(j)* sin(k(j)*b))/ (b * (2*k(j)*h + sin(2*k(j)*h))*(k(j)^2-m(i)^2))...
            *(besselk(l, k(j)*RE)) / (dbesselk(l, k(j)*RE));% Set A_mn
    end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h)); %set B_00 
%%%%%%%%%%%B_0m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
    Rml_prime_RE =  m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) ...
  - besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))...
  /(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
      B(1,j) = Rml_prime_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
 end
%%%%%%%%%%%B_n0%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
    B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));% Set B_n0
end
%%%%%%%%%%%%%%%%%%%%%%%%%%B_nm%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for i=2:N
     for j=2:M
         Rml_prime_RE =  m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) ...
  - besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))...
  /(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
     B(i,j) = Rml_prime_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
      end
  end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * wk*log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for j = 2:M
     Rml_prime_star_RE =  m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) ...
  - besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))...
  /(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
    C(1,j) = Rml_prime_star_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
 end
% %%%%%%% Define Cn0%%%%%%
for i = 2:N
    C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
end
% % %%%%%% Define Cnm%%%%%%
 for i = 2:N
     for j =2:M
 Rml_prime_star_RE =  m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) ...
  - besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))...
  /(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
    C(i,j) = Rml_prime_star_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
     end
 end
%%%%%%% write Matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M,N);
%%% write D00 %%%%%%
D(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2*wk*h + sinh(2*wk*h))) * (besselj(l, wk*RI))/ (dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for j =2:N
    D(1,j) = (cos(k(j)*h) * sin(k(j)*b)) / (k(j)*b * (2*k(j)*h + sin(2*k(j)*h))) * (besseli(l, k(j)*RI) )/ (dbesseli(l, k(j)*RI));

end
%%%%%% write Dm0 %%%%%%%%%
for i = 2:M
    D(i,1) = (2 * cosh(wk*h) * (-1)^(i-1) * wk * sinh(wk*b)) /(b * (2*wk*h + sinh(2*wk*h)) * (wk^2 + m(i)^2))...
        *(besselj(l, wk*RI) )/(dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for i = 2:M
    for j = 2:N
        D(i,j) = (2 * cos(k(j)*h) * (-1)^(i-1) * k(j) * sin(k(j)*b)) /(b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 - m(i)^2))...
        *(besseli(l, k(j)*RI)) / (dbesseli(l, k(j)*RI));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M);  % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
 E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
    Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) ...
  - besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))...
  / (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
 E(1,j) = Rml_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
% 
 end
%%%%%% Define Eno%%%%%%%%%%%%%%
for i = 2:N
    E(i,1) = (4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
for i = 2:N
    for j = 2:M
        Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) ...
  - besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))...
  / (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(i,j) = Rml_prime_RI * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
    end 
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));
% %%%%%% Define F0m%%%%
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) ...
  - besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) ...
  * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(1,j) = Rml_star_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%% Defin Fn0%%%%%%
for i = 2:N
    F(i,1) = (-4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
%%%%%%% Define Fnm %%%%%%%
for i = 2:N
    for j = 2:M
        Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) ...
  - besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) ...
  * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
        F(i,j) = Rml_star_prime_RI * (4 * k(i)* (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
    end
end
% % %%%%%%%%%% Write Matrix W %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 W = zeros(N,Q);
% %%W_{0q}
 for q = 1:Q
     r = exp(-2*chi(q)*b);                          % e^{-2*chi*b}
coth_chib = (1 + r) / (1 - r);                 % coth(chi*b) via decaying exp
 W(1,q) = -(4*chi(q)/cosh(wk*h))*(( wk*sinh(wk*b)*coth_chib - chi(q)*cosh(wk*b) ) / ( wk^2 - chi(q)^2 ) ); %prof.chen
% 
%     % W(1, q) = - (4 * chi(q) / cosh(wk * h)) * (chi(q) * sinh(chi(q) * b) * cosh(wk * b) - wk * sinh(wk * b) * cosh(chi(q) * b)) ...
%     %     / ((chi(q)^2 - wk^2) * sinh(chi(q) * b));
 end
%  % %%% Write W_{nq}
  for n = 2:N
     for q = 1:Q
          r = exp(-2*chi(q)*b);                  % e^{-2*chi*b}
    coth_chib = (1 + r) / (1 - r);         % coth(chi*b), stable
   W(n,q) = -(4*chi(q)/cos(k(n)*h)) * ...
                  ( ( k(n)*sin(k(n)*b)*coth_chib + chi(q)*cos(k(n)*b) ) ...   %prof.chen 
                    / ( k(n)^2 + chi(q)^2 ) );
%         % W(n, q) =  - (4 * chi(q) / cos(k(n) * h)) * (chi(q) * sinh(chi(q) * b) * cos(k(n) * b) + k(n) * sin(k(n) * b) * cosh(chi(q) * b)) ...
%         % / ((chi(q)^2 + k(n)^2) * sinh(chi(q) * b));
      end
  end
 %%%%%%%%%%%%%%%%%%%%%%%Find E1_{q}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 f2 = omega(idx)^2 / g;
 E1  = zeros(Q,1);
S_q = zeros(Q,1);
for q = 1:Q
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% this expression is same as prof. Chen code
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    exp2chi_d = exp(-2 * chi(q) * d);      % e^{-2*chi*d}
    exp2chi_b = exp(-2 * chi(q) * b);      % e^{-2*chi*(h-d)}
    exp2chi_h = exp(-2 * chi(q) * h);      % e^{-2*chi*h}
    exp1chi_d = exp(-chi(q) * d);                 % e^{-chi*d}
    % intermediate ratio (same order as Fortran)
    BB = (f2 - chi(q)) * (exp2chi_b / exp2chi_d);
    BB = BB + (f2 + chi(q)) * (1 + 2 * exp2chi_b);
    BB = BB / (f2 - chi(q) + (f2 + chi(q)) * exp2chi_d);
    BB = BB / (1 + exp2chi_b);
    % final coefficient
    E1(q) = 2 * (BB * exp2chi_d - 1);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% This expression is that one derived directly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 
    %  % N0, D0, H0 exactly as in the math above
    % N0 = (chi(q) - f2) + (chi(q) + f2) * exp2chi_h;
    % D0 = (f2 - chi(q)) + (f2 + chi(q)) * exp2chi_d;
    % H0 = 1 + exp2chi_b;
    % 
    % % final coefficient (no add/subtract step)
    % E1(q) = 2 * N0 / (D0 * H0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This one is the hyperbolic one form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% num = chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h);
% den = (f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*(h - d));  % or: * cosh(chi(q)*b)
% E1(q) = num / den;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% This is the upper region velocity potential at z=0, write
%%%%%%%%%%%%%%%% from prof.chen code 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    CCmc = 1 + (exp2chi_b / exp2chi_d) + 2*exp2chi_b;
 CCmc = CCmc / (1 + exp2chi_b);
 CCmc = BB * (1 + exp2chi_d) - CCmc;
S_q(q)  = CCmc * exp1chi_d;                       % this equals the summand value for this q
% --- direct hyperbolic (exact) ---

% % % C_hyp = ( chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h) ) ...
% % %         / ( (f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*b) );
% % % 
% % % S_hyp(q) = C_hyp * cosh(chi(q)*d) + sinh(chi(q)*h) / cosh(chi(q)*b);   % summand


end


% H(q): diagonal term (exponential form)
H1 = zeros(Q,1);
for q = 1:Q
% dissip = 
    r   = exp(-2*chi(q)*b);                 % e^{-2*chi*b}
    mid   = 4*r/(1 - r^2) - E1(q);           % 2/sinh(2*chi*b) - C_m
    H1(q) = mid + 1i*a1(idx)*chi(q)/omega(idx);            % mid + i*(a1/w)*chi
end
H = diag(H1);
 %%%%%%%%%%matric G%%%%%%%%%%%%%%%%%%
 G = zeros(Q, N);        % Preallocate
 for q = 1:Q
G(q,1) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Z0d ...
       * ( chi(q) / (chi(q)^2 - wk^2) )* ( besselj(l, wk*RI) / dbesselj(l, wk*RI) ); % this one is my calculated 

 end
% %G_qn
for q = 1:Q
    for n = 2:N
         G(q, n) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Znd(n) ...
       * ( chi(q) / (k(n)^2 + chi(q)^2) )* ( besseli(l, k(n)*RI) / dbesseli(l, k(n)*RI) ); % this one is my calculated 

    end
end
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Define the right hand side vector
 U = zeros(2*M + 2*N + Q, 1);  % Full vector, M+N+M+N+Q elements

% % Block 1: U (size M = 4)
for i = 1:M
    if i == 1
        U(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ;  % Z_0^l
    else

        U(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2)*cosh(wk*h))*besselj(l, wk*RE); %Z_m^l
    end
end

% % Block 2: Y (size N = 3)
for j = 1:N
    if j == 1
        U(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2);  % Y_0^l
    else
        U(j + M, 1) = 0;  % Y_n^l
    end
end
% % Block 3: X (size M = 4)
for i = 1:M
    U(i + M + N, 1) = 0;  % X_0^l, X_m^l
end

% Block 4: W (size N = 3)
for j = 1:N
    U(j + M + N + M, 1) = 0;  % W_0^l, W_n^l
end

% % Identity matrices
 I_M = eye(M);
 I_N = eye(N);

 % Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);
ZMQ = zeros(M,Q);
ZNQ = zeros(N,Q);
ZQM = zeros(Q,M);
ZQN = zeros(Q,N);

% Construct the full block matrix S
 S = [ I_M,   -A,     ZMM,  ZMN,  ZMQ;
     -B,       I_N,  -C,   ZNN,  ZNQ;
      ZMM,      ZMN,   I_M, -D,  ZMQ;
     -E,        ZNN,    -F,   I_N, -W;
      ZQM,      ZQN,    ZQM,  -G,   H ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
converged = false;
psi = zeros(Q,1);
T_old = [];                      % <-- track the full solution from prev. iter
% T_old     = zeros(2*M + 2*N + Q, 1);   % previous unknowns (seed)
for iter = 1:max_iter
    % (A) overwrite ONLY Block 5 of U from current psi
    for q = 1:Q
        U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi(q);
    end
  T = S \ U;  % Solve the linear system
 % T = pinv(S) * U;  % Use pseudoinverse for stability

b_vec = T(1:M);                   % Coefficients b^l
a_vec = T(M+1 : M+N);             % Coefficients a^l
c_vec = T(M+N+1 : 2*M+N);         % Coefficients c^l
d_vec = T(2*M+N+1:2*M+2*N);         % Coefficients d^l
e_vec = T(2*M+2*N+1:end);              %  (Q×1)


% (D) update psi for NEXT iteration from CURRENT coefficients
    for q = 1:Q
        integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
                       .* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
                       .* besselj(l, chi(q)*r) .* r;
        psi(q) = integral(integrand, 0, RI, 'AbsTol',1e-8, 'RelTol',1e-6);

    end

% % % drawnow;  % optional: flush output each iter
% === compact per-iteration print, only for T in [5,35] ===
% Tcur = 2*pi/omega(idx);
% if Tcur >= 5 && Tcur <= 35
%     fprintf('idx %3d | iter %2d | ||psi|| = %.3e\n', idx, iter, norm(psi));
%     drawnow;
% end

    % % % % (E) convergence on ALL unknowns (T vs previous T)
   
   % --- per-iteration diagnostics (optional) ---
    if ~isempty(T_old)
    diff_T = max(abs(T - T_old));
    % Tcur = 2*pi/omega(idx);
    % if Tcur >= 5 && Tcur <= 35
    % fprintf('k0 idx %3d | iter %2d: max|ΔT| = %.3e,  ||psi|| = %.3e\n', ...
    %         idx, iter, diff_T, norm(psi));
    % drawnow;  % ← add this to flush each iteration line too
    % end
    if diff_T < tol
        converged = true;
        break
    end
    end
    T_old = T;
end  % <<< end of iter loop

% % ---- ONE summary line per frequency (place it here) ----
% iters_used(idx) = iter;  % how many iterations this frequency needed
% fprintf('idx %3d (T=%6.3f s): iters=%2d\n', idx, 2*pi/omega(idx), iter);
% drawnow;



% -------------- end nonlinear iteration --------------
% % %%%%%%%%% Wave motion%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 term1 = (d_vec(1)*cosh(wk*h)^2) / ((2*wk*h + sinh(2*wk*h)) * dbesselj(l, wk*RI));
sum1 = 0;
for n =2:N
    sum1 = sum1 + d_vec(n)*cos(k(n)*h)^2/ ((2*k(n)*h + sin(2*k(n)*h))* dbesseli(l, k(n)*RI)); 
end 

phiUell = 0;
for q =1:Q
phiUell = phiUell +e_vec(q)* S_q(q);

 end
 eta0(idx) = abs(term1+sum1+phiUell);

   

end


% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %              Plotting Data 
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % 
% % 
% figure(2); clf;
% >>> CHANGES START (3): one figure = two taus + Dr. Chen only



    % plot(2*pi./omega, eta0, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T 
plot(omega.^2 * RE / g, eta0, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T 
hold on;  % Add extra plots on same figure
% % % % 
% % % % % Now plot the 3 CSV experimental points
% scatter(data_chen(:,1), data_chen(:,2), 30, 'r', 'o', 'filled');
%  % scatter(data_liu(:,1),  data_liu(:,2),  30, 'g', 's', 'filled');
      % scatter(data_exp(:,1),  data_exp(:,2),  30, 'b', '^', 'filled');
  xlabel('$T$', 'Interpreter', 'latex');
  ylabel('$|\eta / (iA)|$', 'Interpreter', 'latex');
  title('Wave motion amplitude for $R_E = 138$', 'Interpreter', 'latex');
 legend({'$\tau=0.2$','Model test'}, ...
       'Interpreter','latex','Location','northwest');

  grid on;
  % xlim([5 35]);  % Match the reference plot
  % ylim([0 4.5]);  % Optional: based on expected peak height
xlim([0 4]);  % Match the reference plot
 ylim([0 7]);  % Optional: based on expected peak height
% === plotting section ===


% diary off
elapsedTime = toc;
disp(['Time consuming =  ', num2str(elapsedTime), ' s']);

function out = dbesselk(l, z)
%DBESSELK  Derivative of the modified Bessel function of the second kind
%   out = dbesselk(l, z)
%   Returns d/dz [K_l(z)] using the recurrence formula:
%   K_l'(z) = -0.5 * (K_{l-1}(z) + K_{l+1}(z))

    out = -0.5 * (besselk(l-1, z) + besselk(l+1, z));
end
function out = dbesselj(l, z)
%DBESSELJ  Derivative of the Bessel function of the first kind
%   out = dbesselj(l, z)
%   Returns d/dz [J_l(z)] using the recurrence formula:
%   J_l'(z) = 0.5 * (J_{l-1}(z) - J_{l+1}(z))

    out = 0.5 * (besselj(l-1, z) - besselj(l+1, z));
end
function out = dbesseli(l, z)
%DBESSELI  Derivative of the modified Bessel function of the first kind
%   out = dbesseli(l, z)
%   Returns d/dz [I_l(z)] using the recurrence formula:
%   I_l'(z) = 0.5 * (I_{l-1}(z) + I_{l+1}(z))

    out = 0.5 * (besseli(l-1, z) + besseli(l+1, z));
end





function out = dbesselh(l, z)
%DBESSELH  Derivative of the Hankel function of the first kind
%   out = dbesselh(l, z)
%   Returns d/dz [H_l^{(1)}(z)] using the recurrence formula:
%   H_l^{(1)'}(z) = 0.5 * (H_{l-1}^{(1)}(z) - H_{l+1}^{(1)}(z))

    out = 0.5 * (besselh(l-1, 1, z) - besselh(l+1, 1, z));
end
function x = bessel0j(l,q,opt)
% a row vector of the first q roots of bessel function Jl(x), integer order.
% if opt = 'd', row vector of the first q roots of dJl(x)/dx, integer order.
% if opt is not provided, the default is zeros of Jl(x).
% all roots are positive, except when l=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jl was borrowed from Cleve Moler, 
% but the starting points for both Jl and Jl' can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13. 
%
% David Goodmanson
%
% x = bessel0j(l,q,opt)

k = 1:q;
if nargin==3 && opt=='d'
    beta = (k + l/2 - 3/4)*pi;
    mu = 4*l^2;
    x = beta - (mu+3)./(8*beta) - 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);
    for j=1:8
        xnew = x - besseljd(l,x)./ ...
            (besselj(l,x).*((l^2./x.^2)-1) -besseljd(l,x)./x);
        x = xnew;    
    end
    if l==0
        x(1) = 0;    % correct a small numerical difference from 0
    end
else
    beta = (k + l/2 - 1/4)*pi;
    mu = 4*l^2;
    x = beta - (mu-1)./(8*beta) - 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
    for j=1:8
        xnew = x - besselj(l,x)./besseljd(l,x);
        x = xnew;
    end
end
end

% --- Local helper function for derivative of Bessel function ---
function dJ = besseljd(l, x)
    dJ = 0.5 * (besselj(l - 1, x) - besselj(l + 1, x));
end
function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi)
    % n = 0 mode
    term1 = (d_vec(1) * Z0d * besselj(l, wk*r)) / (dbesselj(l, wk*RI));

    % n >= 2 evanescent sum
    sum1 = 0;
    for nidx = 2:N
        sum1 = sum1 + d_vec(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) ...
                     /( dbesseli(l, k(nidx)*RI));
    end

    % q = 1..Q radial sum
    sum2 = 0;
    for qidx = 1:Q

        sum2 = sum2 + e_vec(qidx) *chi(qidx)* besselj(l, chi(qidx)*r) / dbesselj(l, chi(qidx)*RI);
    end

    v_D_val = term1 + sum1 + sum2;
end

[Javeria]

I was trying to run the code but i get unauthorized to perform this operation however i post my question and code above. for large paramters my code works and i also do validation of this. But when i used the small input parameters then my code give me the numerical instabilities for this v_D non linear function i treated my this one non linear function as iteratively inside my code. How i can solve my this problem that my code also works for this small inpur parameters.

Thanks

[Sam Chak]

Hi @Javeria

Your code is quite lengthy. If it takes more than a minute to run on this platform, you are likely to receive the "unauthorized to perform this operation" message. Have you identified the specific portion of the Bessel function in your code that incurs numerical instabilities? I also notice that part of your code incorporates contributions from @David Goodmanson.

[Javeria]

@Sam Chak my code give me the error for this one portion when i switch to the small input paramaters.

this integral returns me the Inf and NaN values.

(D) update psi for NEXT iteration from CURRENT coefficients
    for q = 1:Q
        integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
                       .* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
                       .* besselj(l, chi(q)*r) .* r;
        psi(q) = integral(integrand, 0, RI, 'AbsTol',1e-8, 'RelTol',1e-6);
    end

and the integral i calculate using this one function,

function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi)
    % n = 0 mode
    term1 = (d_vec(1) * Z0d * besselj(l, wk*r,1)) / (dbesselj(l, wk*RI));

    % n >= 2 evanescent sum
    sum1 = 0;
    for nidx = 2:N
        sum1 = sum1 + d_vec(nidx) * Znd(nidx) * besseli(l, k(nidx)*r,1) ...
                     /( dbesseli(l, k(nidx)*RI));
    end

    % q = 1..Q radial sum
    sum2 = 0;
    for qidx = 1:Q
       
        sum2 = sum2 + e_vec(qidx) *chi(qidx)* besselj(l, chi(qidx)*r,1) / dbesselj(l, chi(qidx)*RI);
    end

    v_D_val = term1 + sum1 + sum2;
end

[Javeria]

i also try the scale one form but same error i get.

[Walter Roberson]

I am not convinced that the problem is due to instability in bessel functions.

        integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
                       .* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
                       .* besselj(l, chi(q)*r) .* r;

after awhile, d_vec and e_vec both grow to roughly +/- 1e172 +/- 1e172i . As a result, v_D returns results in roughly the same range. Multiply two such results together and you get overflow to +/- infinity. The besselj values remain modest.

[Javeria]

@Walter Robersoni check my code just by setting the b1=0 i.e no non linear term and it works for it. Now i am confused that why my code works for large input parameters but not for small input parameters and what might be solution for it.

Thanks

[Sam Chak]

Hi @Javeria

@Walter Roberson noted that both d_vec and e_vec grow to extremely large values and that they are related to the linear system S*T = U. Is the system mathematically guaranteed to have solutions? If so, the very large values in the system are probably causing numerical overflow.

T = S \ U;                      % Solve the linear system
% T = pinv(S) * U;              % Use pseudoinverse for stability

b_vec = T(1:M);                 % Coefficients b^l
a_vec = T(M+1 : M+N);           % Coefficients a^l
c_vec = T(M+N+1 : 2*M+N);       % Coefficients c^l
d_vec = T(2*M+N+1:2*M+2*N);     % Coefficients d^l
e_vec = T(2*M+2*N+1:end);       % (Q×1)

[Javeria]

@Walter Roberson yes the system has solution but how i can fix it for my small input paramters? i attached my this non linear calculation file also. Thanks

[Javeria]

@Walter Roberson@Sam Chak what i have that we need to treat this non linear term iteratively to solve our system of equation;

The principle of iteration is the same. You may go in another way: 1) assuming your coefficient b1=0, you should have everything including v_D(r);

2) using a new a’( r)=a+b|v_D( r)| in place of a to solve all equations to get a new v_D( r); 3) to check the new v_D with previous one, if the difference is larger than your tolerance, continue the step 2).

Answer 1

[Torsten]

Although it's slow, the full Newton method seems to work in this case (except for one value of X_kRE):

clear all;
close all;
clc;
tic;
% diary('iter_log.txt');
% diary on
% fprintf('Run started: %s\n', datestr(now));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%          Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
old = false;
if old
  N = 20;                             % Number of N modes
  M = 20;                             % Number of M modes
  Q = 20;                             % Number of Q modes
  RI = 34.5;                          % Inner radius (m)
  RE = 47.5;                          % Outer radius (m)
  h = 200;                            % Water depth (m) % Interface depth (m)
  d = 38;                             % Interface depth (m) % Draft (m)
  tau = 0.1;                          % porosity ratio
  X_kRE = linspace(0.0025, 0.16,100); %%%%%% this is now k_0
else
  N = 10;                             % Number of N modes
  M = 10;                             % Number of M modes
  Q = 10;                             % Number of Q modes
  RI = 0.4;                           % Inner radius (m)
  RE = 0.8;                           % Outer radius (m)
  h = 1.0;                            % Water depth (m) % Interface depth (m)
  d = 0.3;                            % Interface depth (m) % Draft (m)
  tau = 0.2;                          % porosity ratio
  X_kRE = linspace(0.05, 4.5, 40);    %%%%%% this is now k_0
end
b = h-d;                                             % Distance from cylinder bottom to seabed (m)
g = 9.81;                                            % Gravity (m/s^2)
gamma = 1.0;                                         % discharge coefficient

b1    = (1 - tau) / (2 * gamma * tau^2);   % nonlinear  coeff
tol      = 1e-4;                           % convergence tolerance
max_iter = 100;                            % max  iterations
l = 0;                                     % Azimuthal order
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === Compute Bessel roots and scaled chi values ===
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
roots = bessel0j(l, Q);    % Zeros of J_l(x)
chi = roots ./ RI;         % chi_q^l = roots scaled by radius
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE));  % Preallocate
omega = zeros(size(X_kRE));  %%%% omega preallocate
% iters_used = zeros(length(X_kRE),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for idx = 1:length(X_kRE)
     wk = X_kRE(idx);                                % dimensional wavenumber
        omega(idx) = sqrt(g * wk * tanh(wk * h));    % dispersion relation
        % fprintf('\n--- idx %3d (T=%6.3f s) ---\n', idx, 2*pi/omega(idx));
% drawnow;
%%%%%%%%%%%%% Compute group velocity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Cg(idx) = (g * tanh(wk * h) + g * wk * h * (sech(wk * h))^2) ...
              * omega(idx) / (2 * g * wk * tanh(wk * h));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% =======Compute a1 based on tau
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    a1(idx) = 0.93 * (1 - tau) * Cg(idx);
    % dissip(idx) = a1(idx)/omega(idx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %=============== Find derivative of Z0 at z = -d
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  Z0d = (cosh(wk * h) * wk * sinh(wk * b)) / (2 * wk * h + sinh(2 * wk * h));
  fun_alpha = @(x) omega(idx)^2 + g*x.*tan(x*h);     %dispersion equation for evanescent modes
  opts_lsq = optimset('Display','none');  % ← add this
  %opts_lsq = [];
 for n = 1:N
    if n == 1
        k(n) = -1i * wk;
    else
        x0_left  = (2*n - 3) * pi / (2*h);
        x0_right = (2*n - 1) * pi / (2*h);
        guess = (x0_left + x0_right)/2;
        % Use lsqnonlin for better convergence
         k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right,opts_lsq);
         %%%%%%%%%%%%% Derivative of Z_n at z = -d %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Znd(n) = -k(n) * (cos(k(n)*h) * sin(k(n)*b)) / (2*k(n)*h + sin(2*k(n)*h));

    end
end
%% finding the root \lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;
 end
% %%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  A = zeros(M,N);
   A(1,1) = (cosh(wk*h)*sinh(wk*b))/(wk*b*(2*wk*h+sinh(2*wk*h)))*(besselh(l, wk*RE)) /(dbesselh(l, wk*RE)); % Set A_00
 for j = 2:N
    A(1,j) =(cos(k(j)*h)*sin(k(j)*b) )/ (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))*( besselk(l, k(j)*RE) )...
       /(dbesselk(l, k(j)*RE)); % Set A_0n
 end
for i = 2:M
  A(i,1) =(2 * cosh(wk*h) * (-1)^(i-1) * wk*sinh(wk*b)) / (b * (2*wk*h + sinh(2*wk*h))*(wk^2 + m(i)^2))...
      * (besselh(l, wk*RE))/ (dbesselh(l, wk*RE)); % Set A_m0
end
for i = 2:M
    for j = 2:N
        A(i,j) = (2*cos(k(j)*h)*(-1)^(i-1) *k(j)* sin(k(j)*b))/ (b * (2*k(j)*h + sin(2*k(j)*h))*(k(j)^2-m(i)^2))...
            *(besselk(l, k(j)*RE)) / (dbesselk(l, k(j)*RE));% Set A_mn
    end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h)); %set B_00
%%%%%%%%%%%B_0m%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
    Rml_prime_RE =  m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) ...
  - besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))...
  /(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
      B(1,j) = Rml_prime_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
 end
%%%%%%%%%%%B_n0%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N
    B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));% Set B_n0
end
%%%%%%%%%%%%%%%%%%%%%%%%%%B_nm%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 for i=2:N
     for j=2:M
         Rml_prime_RE =  m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RE) ...
  - besseli(l, m(j)*RI) * dbesselk(l, m(j)*RE))...
  /(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
     B(i,j) = Rml_prime_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
      end
  end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * wk*log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for j = 2:M
     Rml_prime_star_RE =  m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) ...
  - besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))...
  /(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
    C(1,j) = Rml_prime_star_RE * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
 end
% %%%%%%% Define Cn0%%%%%%
for i = 2:N
    C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
end
% % %%%%%% Define Cnm%%%%%%
 for i = 2:N
     for j =2:M
 Rml_prime_star_RE =  m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RE) ...
  - besselk(l, m(j)*RE) * dbesseli(l, m(j)*RE))...
  /(besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
    C(i,j) = Rml_prime_star_RE * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
     end
 end
%%%%%%% write Matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M,N);
%%% write D00 %%%%%%
D(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2*wk*h + sinh(2*wk*h))) * (besselj(l, wk*RI))/ (dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for j =2:N
    D(1,j) = (cos(k(j)*h) * sin(k(j)*b)) / (k(j)*b * (2*k(j)*h + sin(2*k(j)*h))) * (besseli(l, k(j)*RI) )/ (dbesseli(l, k(j)*RI));
end
%%%%%% write Dm0 %%%%%%%%%
for i = 2:M
    D(i,1) = (2 * cosh(wk*h) * (-1)^(i-1) * wk * sinh(wk*b)) /(b * (2*wk*h + sinh(2*wk*h)) * (wk^2 + m(i)^2))...
        *(besselj(l, wk*RI) )/(dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for i = 2:M
    for j = 2:N
        D(i,j) = (2 * cos(k(j)*h) * (-1)^(i-1) * k(j) * sin(k(j)*b)) /(b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 - m(i)^2))...
        *(besseli(l, k(j)*RI)) / (dbesseli(l, k(j)*RI));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M);  % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
 E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 2:M
    Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) ...
  - besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))...
  / (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
 E(1,j) = Rml_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
%
 end
%%%%%% Define Eno%%%%%%%%%%%%%%
for i = 2:N
    E(i,1) = (4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
for i = 2:N
    for j = 2:M
        Rml_prime_RI = m(j)*(besselk(l, m(j)*RI) * dbesseli(l, m(j)*RI) ...
  - besseli(l, m(j)*RI) * dbesselk(l, m(j)*RI))...
  / (besselk(l, m(j)*RI) * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI));
E(i,j) = Rml_prime_RI * (4 * k(i) * (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));
% %%%%%% Define F0m%%%%
for j = 2:M
Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) ...
  - besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) ...
  * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
F(1,j) = Rml_star_prime_RI * (4 * wk * (-1)^(j-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + m(j)^2));
end
%%%%% Defin Fn0%%%%%%
for i = 2:N
    F(i,1) = (-4 * sin(k(i)*b)) / (RI *k(i)* log(RE/RI) * cos(k(i)*h));
end
%%%%%%% Define Fnm %%%%%%%
for i = 2:N
    for j = 2:M
        Rml_star_prime_RI = m(j)*(besseli(l, m(j)*RE) * dbesselk(l, m(j)*RI) ...
  - besselk(l, m(j)*RE) * dbesseli(l, m(j)*RI))/((besselk(l, m(j)*RI) ...
  * besseli(l, m(j)*RE) - besselk(l, m(j)*RE) * besseli(l, m(j)*RI)));
        F(i,j) = Rml_star_prime_RI * (4 * k(i)* (-1)^(j-1) * sin(k(i)*b)) / (cos(k(i)*h) * (k(i)^2 - m(j)^2));
    end
end
% % %%%%%%%%%% Write Matrix W %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 W = zeros(N,Q);
% %%W_{0q}
 for q = 1:Q
     r = exp(-2*chi(q)*b);                          % e^{-2*chi*b}
coth_chib = (1 + r) / (1 - r);                 % coth(chi*b) via decaying exp
 W(1,q) = -(4*chi(q)/cosh(wk*h))*(( wk*sinh(wk*b)*coth_chib - chi(q)*cosh(wk*b) ) / ( wk^2 - chi(q)^2 ) ); %prof.chen
%
%     % W(1, q) = - (4 * chi(q) / cosh(wk * h)) * (chi(q) * sinh(chi(q) * b) * cosh(wk * b) - wk * sinh(wk * b) * cosh(chi(q) * b)) ...
%     %     / ((chi(q)^2 - wk^2) * sinh(chi(q) * b));
 end
%  % %%% Write W_{nq}
  for n = 2:N
     for q = 1:Q
          r = exp(-2*chi(q)*b);                  % e^{-2*chi*b}
    coth_chib = (1 + r) / (1 - r);         % coth(chi*b), stable
   W(n,q) = -(4*chi(q)/cos(k(n)*h)) * ...
                  ( ( k(n)*sin(k(n)*b)*coth_chib + chi(q)*cos(k(n)*b) ) ...   %prof.chen
                    / ( k(n)^2 + chi(q)^2 ) );
%         % W(n, q) =  - (4 * chi(q) / cos(k(n) * h)) * (chi(q) * sinh(chi(q) * b) * cos(k(n) * b) + k(n) * sin(k(n) * b) * cosh(chi(q) * b)) ...
%         % / ((chi(q)^2 + k(n)^2) * sinh(chi(q) * b));
      end
  end
 %%%%%%%%%%%%%%%%%%%%%%%Find E1_{q}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 f2 = omega(idx)^2 / g;
 E1  = zeros(Q,1);
S_q = zeros(Q,1);
for q = 1:Q
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% this expression is same as prof. Chen code
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    exp2chi_d = exp(-2 * chi(q) * d);      % e^{-2*chi*d}
    exp2chi_b = exp(-2 * chi(q) * b);      % e^{-2*chi*(h-d)}
    exp2chi_h = exp(-2 * chi(q) * h);      % e^{-2*chi*h}
    exp1chi_d = exp(-chi(q) * d);                 % e^{-chi*d}
    % intermediate ratio (same order as Fortran)
    BB = (f2 - chi(q)) * (exp2chi_b / exp2chi_d);
    BB = BB + (f2 + chi(q)) * (1 + 2 * exp2chi_b);
    BB = BB / (f2 - chi(q) + (f2 + chi(q)) * exp2chi_d);
    BB = BB / (1 + exp2chi_b);
    % final coefficient
    E1(q) = 2 * (BB * exp2chi_d - 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% This expression is that one derived directly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    %  % N0, D0, H0 exactly as in the math above
    % N0 = (chi(q) - f2) + (chi(q) + f2) * exp2chi_h;
    % D0 = (f2 - chi(q)) + (f2 + chi(q)) * exp2chi_d;
    % H0 = 1 + exp2chi_b;
    %
    % % final coefficient (no add/subtract step)
    % E1(q) = 2 * N0 / (D0 * H0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This one is the hyperbolic one form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% num = chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h);
% den = (f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*(h - d));  % or: * cosh(chi(q)*b)
% E1(q) = num / den;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% This is the upper region velocity potential at z=0, write
%%%%%%%%%%%%%%%% from prof.chen code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    CCmc = 1 + (exp2chi_b / exp2chi_d) + 2*exp2chi_b;
 CCmc = CCmc / (1 + exp2chi_b);
 CCmc = BB * (1 + exp2chi_d) - CCmc;
S_q(q)  = CCmc * exp1chi_d;                       % this equals the summand value for this q
% --- direct hyperbolic (exact) ---
% % % C_hyp = ( chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h) ) ...
% % %         / ( (f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d)) * cosh(chi(q)*b) );
% % %
% % % S_hyp(q) = C_hyp * cosh(chi(q)*d) + sinh(chi(q)*h) / cosh(chi(q)*b);   % summand
end
% H(q): diagonal term (exponential form)
H1 = zeros(Q,1);
for q = 1:Q
% dissip =
    r   = exp(-2*chi(q)*b);                 % e^{-2*chi*b}
    mid   = 4*r/(1 - r^2) - E1(q);           % 2/sinh(2*chi*b) - C_m
    H1(q) = mid + 1i*a1(idx)*chi(q)/omega(idx);            % mid + i*(a1/w)*chi
end
H = diag(H1);
 %%%%%%%%%%matric G%%%%%%%%%%%%%%%%%%
 G = zeros(Q, N);        % Preallocate
 for q = 1:Q
G(q,1) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Z0d ...
       * ( chi(q) / (chi(q)^2 - wk^2) )* ( besselj(l, wk*RI) / dbesselj(l, wk*RI) ); % this one is my calculated
 end
% %G_qn
for q = 1:Q
    for n = 2:N
         G(q, n) = 1i * 2*(a1(idx)/(omega(idx)*RI))* Znd(n) ...
       * ( chi(q) / (k(n)^2 + chi(q)^2) )* ( besseli(l, k(n)*RI) / dbesseli(l, k(n)*RI) ); % this one is my calculated
    end
end
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %Define the right hand side vector
 U = zeros(2*M + 2*N + Q, 1);  % Full vector, M+N+M+N+Q elements
% % Block 1: U (size M = 4)
for i = 1:M
    if i == 1
        U(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ;  % Z_0^l
    else
        U(i, 1) = (2 *wk*(-1)^(i-1)*sinh(wk*b))/(b *(wk^2 +m(i)^2)*cosh(wk*h))*besselj(l, wk*RE); %Z_m^l
    end
end
% % Block 2: Y (size N = 3)
for j = 1:N
    if j == 1
        U(j + M, 1) = -dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2);  % Y_0^l
    else
        U(j + M, 1) = 0;  % Y_n^l
    end
end
% % Block 3: X (size M = 4)
for i = 1:M
    U(i + M + N, 1) = 0;  % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
    U(j + M + N + M, 1) = 0;  % W_0^l, W_n^l
end

for q = 1:Q
   U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI)));
end

% % Identity matrices
 I_M = eye(M);
 I_N = eye(N);
 % Zero matrices
ZMM = zeros(M, M);
ZMN = zeros(M, N);
ZNM = zeros(N, M);
ZNN = zeros(N, N);
ZMQ = zeros(M,Q);
ZNQ = zeros(N,Q);
ZQM = zeros(Q,M);
ZQN = zeros(Q,N);
% Construct the full block matrix S
 S = [ I_M,   -A,     ZMM,  ZMN,  ZMQ;
     -B,       I_N,  -C,   ZNN,  ZNQ;
      ZMM,      ZMN,   I_M, -D,  ZMQ;
     -E,        ZNN,    -F,   I_N, -W;
      ZQM,      ZQN,    ZQM,  -G,   H ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
solution_method = 0;
if solution_method == 0

  if idx==1
    psi0 = zeros(Q,1);
  else
    if exitflag == 1
      psi0 = psi;
    else
      psi0 = zeros(Q,1);
    end
  end
  f = @(psi)fun(psi,M,N,Q,S,U,Z0d,wk,RI,l,Znd,k,chi);
  [psi,~,exitflag] = fsolve(f,psi0,optimset('Display','none'));
  [~,d_vec,e_vec] = f(psi);

elseif solution_method == 1

  converged = false;
  psi = zeros(Q,1);
  T_old = [];                      % <-- track the full solution from prev. iter
  % T_old     = zeros(2*M + 2*N + Q, 1);   % previous unknowns (seed)
  for iter = 1:max_iter
    % (A) overwrite ONLY Block 5 of U from current psi
    for q = 1:Q
        U(2*M + 2*N + q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi(q);
    end
    T = S \ U;  % Solve the linear system
   % T = pinv(S) * U;  % Use pseudoinverse for stability
    b_vec = T(1:M);                   % Coefficients b^l
    a_vec = T(M+1 : M+N);             % Coefficients a^l
    c_vec = T(M+N+1 : 2*M+N);         % Coefficients c^l
    d_vec = T(2*M+N+1:2*M+2*N);         % Coefficients d^l
    e_vec = T(2*M+2*N+1:end);              %  (Q×1)
    % (D) update psi for NEXT iteration from CURRENT coefficients
    for q = 1:Q
        integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
                       .* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
                       .* besselj(l, chi(q)*r) .* r;
        psi(q) = integral(integrand, 0, RI, 'AbsTol',1e-8, 'RelTol',1e-6);
    end
    % % % drawnow;  % optional: flush output each iter
    % === compact per-iteration print, only for T in [5,35] ===
    % Tcur = 2*pi/omega(idx);
    % if Tcur >= 5 && Tcur <= 35
    %     fprintf('idx %3d | iter %2d | ||psi|| = %.3e\n', idx, iter, norm(psi));
    %     drawnow;
    % end
    % % % % (E) convergence on ALL unknowns (T vs previous T)

   % --- per-iteration diagnostics (optional) ---
    if ~isempty(T_old)
      diff_T = max(abs(T - T_old));
      % Tcur = 2*pi/omega(idx);
      % if Tcur >= 5 && Tcur <= 35
      % fprintf('k0 idx %3d | iter %2d: max|ΔT| = %.3e,  ||psi|| = %.3e\n', ...
      %         idx, iter, diff_T, norm(psi));
      % drawnow;  % ← add this to flush each iteration line too
      % end
      if diff_T < tol
        converged = true;
        break
      end
    end
    T_old = T;
  end  % <<< end of iter loop
end
if solution_method==0
  exitflag
elseif solution_method==1
  converged
end
%T
% % ---- ONE summary line per frequency (place it here) ----
% iters_used(idx) = iter;  % how many iterations this frequency needed
% fprintf('idx %3d (T=%6.3f s): iters=%2d\n', idx, 2*pi/omega(idx), iter);
% drawnow;
% -------------- end nonlinear iteration --------------
% % %%%%%%%%% Wave motion%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 term1 = (d_vec(1)*cosh(wk*h)^2) / ((2*wk*h + sinh(2*wk*h)) * dbesselj(l, wk*RI));
sum1 = 0;
for n =2:N
    sum1 = sum1 + d_vec(n)*cos(k(n)*h)^2/ ((2*k(n)*h + sin(2*k(n)*h))* dbesseli(l, k(n)*RI));
end
phiUell = 0;
for q =1:Q
phiUell = phiUell +e_vec(q)* S_q(q);
 end
 eta0(idx) = abs(term1+sum1+phiUell);

end
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %              Plotting Data
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% figure(2); clf;
% >>> CHANGES START (3): one figure = two taus + Dr. Chen only
%if old
  %plot(2*pi./omega, eta0, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T
%else
  plot(omega.^2 * RE / g, eta0, 'k', 'LineWidth', 1.5); %%% this is T which is T= 2*pi/omega plotting against T
%end
hold on;  % Add extra plots on same figure
% % % %
% % % % % Now plot the 3 CSV experimental points
% scatter(data_chen(:,1), data_chen(:,2), 30, 'r', 'o', 'filled');
%  % scatter(data_liu(:,1),  data_liu(:,2),  30, 'g', 's', 'filled');
      % scatter(data_exp(:,1),  data_exp(:,2),  30, 'b', '^', 'filled');
  xlabel('$T$', 'Interpreter', 'latex');
  ylabel('$|\eta / (iA)|$', 'Interpreter', 'latex');
  title('Wave motion amplitude for $R_E = 138$', 'Interpreter', 'latex');
 legend({'$\tau=0.2$','Model test'}, ...
       'Interpreter','latex','Location','northwest');
  grid on;
  % xlim([5 35]);  % Match the reference plot
  % ylim([0 4.5]);  % Optional: based on expected peak height
%xlim([0 4]);  % Match the reference plot
% ylim([0 7]);  % Optional: based on expected peak height
% === plotting section ===
% diary off
elapsedTime = toc;
disp(['Time consuming =  ', num2str(elapsedTime), ' s']);

function [res,d_vec,e_vec] = fun(psi,M,N,Q,S,U,Z0d,wk,RI,l,Znd,k,chi)
  for q = 1:Q
    U(2*M + 2*N + q) = U(2*M + 2*N + q) * psi(q);
  end
  T = S\U;
  d_vec = T(2*M+N+1:2*M+2*N);         % Coefficients d^l
  e_vec = T(2*M+2*N+1:end);           %  (Q×1)
  for q = 1:Q
    integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
                     .* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
                     .* besselj(l, chi(q)*r) .* r;
    res(q) = psi(q) - integral(integrand, 0, RI, 'AbsTol',1e-8, 'RelTol',1e-6);
  end
end
  
function out = dbesselk(l, z)
%DBESSELK  Derivative of the modified Bessel function of the second kind
%   out = dbesselk(l, z)
%   Returns d/dz [K_l(z)] using the recurrence formula:
%   K_l'(z) = -0.5 * (K_{l-1}(z) + K_{l+1}(z))
    out = -0.5 * (besselk(l-1, z) + besselk(l+1, z));
end
function out = dbesselj(l, z)
%DBESSELJ  Derivative of the Bessel function of the first kind
%   out = dbesselj(l, z)
%   Returns d/dz [J_l(z)] using the recurrence formula:
%   J_l'(z) = 0.5 * (J_{l-1}(z) - J_{l+1}(z))
    out = 0.5 * (besselj(l-1, z) - besselj(l+1, z));
end
function out = dbesseli(l, z)
%DBESSELI  Derivative of the modified Bessel function of the first kind
%   out = dbesseli(l, z)
%   Returns d/dz [I_l(z)] using the recurrence formula:
%   I_l'(z) = 0.5 * (I_{l-1}(z) + I_{l+1}(z))
    out = 0.5 * (besseli(l-1, z) + besseli(l+1, z));
end
function out = dbesselh(l, z)
%DBESSELH  Derivative of the Hankel function of the first kind
%   out = dbesselh(l, z)
%   Returns d/dz [H_l^{(1)}(z)] using the recurrence formula:
%   H_l^{(1)'}(z) = 0.5 * (H_{l-1}^{(1)}(z) - H_{l+1}^{(1)}(z))
    out = 0.5 * (besselh(l-1, 1, z) - besselh(l+1, 1, z));
end
function x = bessel0j(l,q,opt)
% a row vector of the first q roots of bessel function Jl(x), integer order.
% if opt = 'd', row vector of the first q roots of dJl(x)/dx, integer order.
% if opt is not provided, the default is zeros of Jl(x).
% all roots are positive, except when l=0,
% x=0 is included as a root of dJ0(x)/dx (standard convention).
%
% starting point for for zeros of Jl was borrowed from Cleve Moler,
% but the starting points for both Jl and Jl' can be found in
% Abramowitz and Stegun 9.5.12, 9.5.13.
%
% David Goodmanson
%
% x = bessel0j(l,q,opt)
k = 1:q;
if nargin==3 && opt=='d'
    beta = (k + l/2 - 3/4)*pi;
    mu = 4*l^2;
    x = beta - (mu+3)./(8*beta) - 4*(7*mu^2+82*mu-9)./(3*(8*beta).^3);
    for j=1:8
        xnew = x - besseljd(l,x)./ ...
            (besselj(l,x).*((l^2./x.^2)-1) -besseljd(l,x)./x);
        x = xnew;
    end
    if l==0
        x(1) = 0;    % correct a small numerical difference from 0
    end
else
    beta = (k + l/2 - 1/4)*pi;
    mu = 4*l^2;
    x = beta - (mu-1)./(8*beta) - 4*(mu-1)*(7*mu-31)./(3*(8*beta).^3);
    for j=1:8
        xnew = x - besselj(l,x)./besseljd(l,x);
        x = xnew;
    end
end
end
% --- Local helper function for derivative of Bessel function ---
function dJ = besseljd(l, x)
    dJ = 0.5 * (besselj(l - 1, x) - besselj(l + 1, x));
end
function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi)
    % n = 0 mode
    term1 = (d_vec(1) * Z0d * besselj(l, wk*r)) / (dbesselj(l, wk*RI));
    % n >= 2 evanescent sum
    sum1 = 0;
    for nidx = 2:N
        sum1 = sum1 + d_vec(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) ...
                     /( dbesseli(l, k(nidx)*RI));
    end
    % q = 1..Q radial sum
    sum2 = 0;
    for qidx = 1:Q
        sum2 = sum2 + e_vec(qidx) *chi(qidx)* besselj(l, chi(qidx)*r) / dbesselj(l, chi(qidx)*RI);
    end
    v_D_val = term1 + sum1 + sum2;
end

[Javeria]

@Torsten Thanks Torsten , As you mentioned that (Although it's slow, the full Newton method seems to work in this case (except for one value of X_kRE)) , our solution must converged for each value of X_kRE did you mean that i need to modify this range??

[Torsten]

did you mean that i need to modify this range??

No. If you run the code, you will see what I mean.

I'd just remove that point in the plot. Or maybe you find a way to make the code converge for this 1 value out of the 40 you prescribed, too.

[Javeria]

@Torsten by running your code i get this plot.

[Javeria]

@Torsten what is the main cause of error in my code that i get diverge solution.

[Javeria]

@Torsten@Walter Roberson@Sam Chak Is this numerical instabilities is due to the bessel functions? I read about that then we used the scaled version of bessel functions. I just make use of it like besseli(n,z,scale) but this way did not solved my problem.

[Torsten]

Your iteration is some kind of fixed-point iteration, Newton's method uses derivative information to get the roots of the nonlinear system. So, Newton's method uses more information about the functions in its iteration process and thus usually has a larger region of convergence than the fixed-point iteration. Further, "fsolve" (the nonlinear solver I used) uses a damped Newton method that is even more sophisticated than the classical Newton method. I think it will lead you nowhere to speculate about a particular reason for divergence in your special case.

And you shouldn't interprete too much into the "second peak" in the result's graph. If you look at the exitflag from the nonlinear solver, it changes here from 1 to 3 (I think) in the point before the second peak (and has a lot of trouble converging in the points nearby as you can see from the long iteration times for these points). So I guess the result before the second peak (the T ~ 0.8 result) is not trustworthy, and at the moment I don't know how to make the solver converge for this particular value of X_kRE.

Because I start the iteration in point i from a trustworthy result in step (i-1), I tried with X_kRE = linspace(0.05, 4.5, 80); instead of X_kRE = linspace(0.05, 4.5, 40); to be nearer to the new solution when restarting, but this even gives exitflag = -3 for two points in the critical region. For the meaning of the exitflags, see

https://uk.mathworks.com/help/optim/ug/fsolve.html

[Javeria]

@Torsten I am curious that by using the large input parameters and use this one

X_kRE = linspace(0.0025, 0.16,100);

i get the converged solution i write this code using these input parameters(deep water condition). But when i switch to the small input parameters(intermediate water depth) then i get this error.

[Torsten]

You get a converged solution for both cases - with the exception of one input point - if you use a sophisticated root finder as in the code above.

I tried to explain that fixed-point methods usually have a smaller radius of convergence around the solution than Newton-like methods.

Since both methods (your method and Newton's method) start from the zero solution and your fixed-point method fails while Newton's method succeeds, this is the most probable reason for the different behaviour.

The physics for the intermediate water depth case seems to be numerically more challenging than the deep water case and needs a more sophisticated solver for the nonlinear contribution. Maybe you find a reason in the physics of the two problems why this is the case. If you look at the problem only from the numerical standpoint, it's too complex to identify the reason for failure with your solution method (at least for me).

Numerics is trial and error: if you know that your equations are correct and your solution method fails, you usually try another solution method. If this works, you don't ask why the first failed because it usually cannot be answered easily.

By the way: Newton's method also works for the deep water conditions (although it's again slower than fixed-point iteration which also succeeds in this case).

[Javeria]

@Torsten Thanks for you kind suggestions i will look it.

[Sam Chak]

Hi @Javeria

I couldn't help because as I haven't followed your previous discussions. However, if your studies focus on wave motion, your aim should be to select the most effective and efficient solver. If your research involves successfully applying a novel numerical method to solve wave motion, you need not abandon the proposed method simply because it is challenging to solve the system at intermediate water depths.

In fact, Newtonian mechanics, which assumes that an object's mass remains constant regardless of its velocity (𝑣), also fails at speeds approaching the speed of light (𝑐). In reality, mass increases as an object approaches the speed of light, approaching infinity as 𝑣 nears 𝑐.

Moreover, your fixed-point iteration method outperforms Newton's method in the deep water case. Therefore, you may consider a hybrid solver approach that mimics Continuously Variable Transmission (CVT) smoothness in automatically selecting Newton's method for intermediate water depths and the fixed-point iteration method for the deep water case.

[Javeria]

@Sam Chak Okay thanks Sam, however if we choose chebeshev guass quadrature for solving this problem with out changing the other things just to deal with this non-linear term.

[Torsten]

@Javeria

I don't think that changing the code to compute the integral

        integrand = @(r) abs(v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi)) ...
                       .* v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi) ...
                       .* besselj(l, chi(q)*r) .* r;
        psi(q) = integral(integrand, 0, RI, 'AbsTol',1e-8, 'RelTol',1e-6);

can help to solve the present problem.

What you have to do is to solve a system of equations

S*T = U(T)

for T.

What you do in your code is to perform a fixed-point iteration

T_(n+1) = S\U(T_n)

and hope that T_n will converge to some vector T.

The equivalent formulation would be to use Newton's method to solve

S*T - U(T) = 0

for T (I do something slightly different since I don't solve for T, but for psi).

[Javeria]

@Torsten Thsnks but at the end we have the converged solution for each X_kRE(this is actually k0 in my code). I attached some refrences.

[Torsten]

Out of interest: how did you change your code to get the iteration converge ?

[Javeria]

@Torsten then how i can prevent from this numerically instabilties for small input parameters? I want converged solution for each k0 not for someone.

Thanks

[David Goodmanson]

Hi Sam,

The concept of equivalent relativistic mass is ok for colloquial arguments, but physics has kind of moved away from that concept. I think that's mostly because of the inclusion of general relativity into the mainstream, where people think about spacetime and other funcamental quantities like (E,p) that transform relativistically. If m is the invariant rest mass, then

E = mc^2/sqrt(1-v^2/c^2)          
p =   mv/sqrt(1-v^2/c^2)              

E^2 = p^2c^2 + m^2c^4
   p = 0   ==>  E = mc^2          % object at rest   

It's algebraically correct but a bit dubious to conclude from the first and last lines that

E/c^2 = m/sqrt(1-v^2/c^2)  = m_equiv.

m_equiv is a quantity that kind of sits on its own. Without it, the speed-up argument is the same but different: To take a massive (rest mass) object up to v=c in a given inertial frame you have to pump in an unllimited amount of energy.

[Torsten]

@Javeria

Thsnks but at the end we have the converged solution for each X_kRE(this is actually k0 in my code).

Ok, I misunderstood this comment from your side. I interpreted it that you already solved your problem with the help of the attached article.

Which part of your code could be affected by the article ?

[Javeria]

@Torsten they treated only the non linear integral with chebeshev guass quadrature and then used the damped newton method.

[Torsten]

I used the modified code to compute eta0 as a function of X_kRE with a finer resolution around the peak. I got the following result:

A = readmatrix('data.txt');
plot(A(:,1),A(:,2))
xlabel('$T$', 'Interpreter', 'latex');
ylabel('$|\eta / (iA)|$', 'Interpreter', 'latex');
title('Wave motion amplitude for $R_E = 138$', 'Interpreter', 'latex');
legend({'$\tau=0.2$'}, ...
       'Interpreter','latex','Location','northwest');
grid on

●

[Javeria]

@Torsten thanks could you please send the full script.

[Torsten]

I left your computational part unchanged, but had to add some code in order to make it work according to my ideas.

In the additionally attached "program1.m", you can also use your fixed-point iteration if you set "solution_method = 0". This option will not work in "program.m". So use "program1.m" instead of "program.m" since it is an extension of "program.m".

If you have questions, you can of course come back and ask.