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).

[David Goodmanson]

Hi Sam, thanks for mentioning that I wrote a piece of the code, which I posted to Javeria on a previous thread. Javeria, if you see this, you can use the bessel0j code as you see fit but the nitpicky part of my nature wishes you had let the variabIe n (the bessel function order) alone, instead of changing it to l (small L). It's just a detail but in code generally, there are too many fonts where small L and capital i and the number 1 can be confused.

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.

[Javeria]

@Torsten Thanks

[Javeria]

@Torsten hello i checked your code can you a little explain that which method you used in the program1 to solve the non linear integral and what modification you did for it.

thanks

[Javeria]

@Torsten in the above code programe1 you write

old = true; %%%%%%%%%% change whn need 
  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)
    d = 38;                             % Interface depth (m) % Draft (m)
    tau = 0.1;                          % porosity ratio
    X_kRE_min = 0.0025;
    X_kRE_max = 0.16;                   %%%%%% this is now k_0 ;
    nXsteps = 100;
    stepX = (X_kRE_max-X_kRE_min)/...
             nXsteps;
    X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
    stepX_min = stepX/2^6;
  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_min = 2.0653;%0.05;%2.0653;
    X_kRE_max = 2.5914;%4.5;%2.5914;            %%%%%% this is now k_0 ;
    nXsteps = 160;%40;%160;
    stepX = (X_kRE_max-X_kRE_min)/...
             nXsteps;
    X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
    stepX_min = stepX/2^15;

what i get that this code is running for the large scale parameters, i change just this one line to false instead of true i get the follwoing error

Warning: Infinite or Not-a-Number value encountered.

> In quadgk>vadapt (line 281)

In quadgk (line 184)

In program1>get_solution (line 592)

In program1>main (line 78)

In program1 (line 1)

Unrecognized function or variable 'T'.

Error in program1>main (line 164)

AA = [T.',Eta0.'];

Error in program1 (line 1)

main()

    old = true; %%%%%%%%%% change whn need 

[Torsten]

If you change "ode=true" to "old=false" and leave "solution_method = 0", you try to run your fixed-point iteration method together with these parameters:

    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
    

This will not work and that's why you came here to ask for a solution.

Set "solution_method" to 1 or 2 to make it work.

Concerning the integration method, I use "quadgk", a Gauss-Kronos quadrature rule, instead of "integral". In some cases, this worked better than "integral". If you want to return to "integral", comment out the lines where "quadgk" is called and replace them by the "integral" calls (that are still present at these code positions).

I additionally added a simple trapezoidal rule to evaluate the integral which can be chosen by using "method_integral = 1" and an appropriate number of discretization points "nr".

If you have further questions, please include the complete parameter setting part:

  solution_method = 0;                  % 0: Fixed-point iteration
                                        % 1: Newton method (psi-based)
                                        % 2: Newton method (minimal-T-based)
                                        % 3: Newton method (medium-T-based) (too slow, not recommended)
                                        % 4: Newton method (full-T-based) (too slow, not recommended)
  method_integral = 0;                  % 0: Use integral
                                        % 1: Use trapezoidal rule
  nr = 300;                             % if method_integral = 1: number of discretization points in r-direction

  % Used if method_integral = 0
  RelTol_integral = 1e-6;               % Relative error tolerance in integral solving (default value: 1e-6)
  AbsTol_integral = 1e-10;              % Absolute error tolerance in integral solving (default value: 1e-10)

  % Used if solution_method = 0
  urelax = 0.01;                        % Underrelaxation factor
  tol = 1e-4;                           % convergence tolerance
  max_iter = 500;                       % max  iterations

  % Used if solution_method ~= 0
  TolX_nleq = 1e-6;                     % Relative error tolerance in integral solving (default value: 1e-6)
  TolFun_nleq = 1e-6;                   % Absolute error tolerance in integral solving (default value: 1e-6)
  AutoScaling_nleq = 'off';             % default value: 'off'
  ComplexEqn_nleq = 'off';              % default value: 'off'

  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)
    d = 38;                             % Interface depth (m) % Draft (m)
    tau = 0.1;                          % porosity ratio
    X_kRE_min = 0.0025;
    X_kRE_max = 0.16;                   %%%%%% this is now k_0 ;
    nXsteps = 100;
    stepX = (X_kRE_max-X_kRE_min)/...
             nXsteps;
    X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
    stepX_min = stepX/2^6;
  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_min = 2.0653;%0.05;%2.0653;
    X_kRE_max = 2.5914;%4.5;%2.5914;            %%%%%% this is now k_0 ;
    nXsteps = 160;%40;%160;
    stepX = (X_kRE_max-X_kRE_min)/...
             nXsteps;
    X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
    stepX_min = stepX/2^15;
  end
  

Usually, it should not be necessary to change parts of the code below these lines and so it will help to identify the problem.

[Javeria]

@Torsten i make the changes according to you as you mentioned above and check the code program1 but i am facing these error

>> program1

Unrecognized function or variable 'T'.

Error in program1>main (line 164)

AA = [T.',Eta0.'];

Error in program1 (line 1)

main()

also i run your code programm only and it takes more than 17 hours but i did not get any result the code is still in running.

[Torsten]

I developed the code under Octave, and the exitflags of "fsolve" under octave are different. So please only use the below code with solver_method = 11 or solver_method = 12 as these options invoked code that is "fsolve" - independent. But I strongly recommend to make the "fsolve" option work under MATLAB, too (the code below works without problems under Octave with solver_method = 21 and solver_method = 22):

https://octave.org/

Under MATLAB answers, the code below needed 32 seconds for 6 values of the X_kRE array. I requested 161 - thus the approximate runtime should be 32*161/6 seconds (maybe some more, because more iterations are needed for bigger values of X_kRE). On my Laptop, it took 3016 sec under octave. But you can see the progress on the terminal.

Please run the below code and report if it works.

After this, remove the lines

    nXsteps = 5;
    X_kRE=X_kRE(1:6);

to run the full second case with 160 instead of 6 output points in the critical region. How long does it take on your computer ?

main()

iter = 8

Starting computations ...
 

iter = 3

x = 2.0686

dx = 0.0033

eta0 = 1.7356

iter = 3

x = 2.0719

dx = 0.0033

eta0 = 1.7526

iter = 3

x = 2.0752

dx = 0.0033

eta0 = 1.7701

iter = 3

x = 2.0785

dx = 0.0033

eta0 = 1.7881

iter = 3

x = 2.0817

dx = 0.0033

eta0 = 1.8065

ans = 1×2

1 6

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ans = 1×2

1 6

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Elapsed time is 31.914075 seconds.

function main

  tic
  %warning('off')
  solution_method = 2;                  % -4: Fixed-point iteration (full-T-based) (unstable, not recommended)
                                        % -3: Fixed-point iteration (medium-T-based) (unstable, not recommended)
                                        % -2: Fixed-point iteration (minimal-T-based) (unstable, not recommended)
                                        % -1: Fixed-point iteration (psi-based) (unstable, not recommended)
                                        % 0: Fixed-point iteration (Ali) (unstable, not recommended)
                                        % 1: Newton method (psi-based)
                                        % 2: Newton method (minimal-T-based)
                                        % 3: Newton method (medium-T-based) (too slow, not recommended)
                                        % 4: Newton method (full-T-based) (too slow, not recommended)
  solver_method = 12;                   % if solution_method > 0:
                                        % 11: Use self-programmed Newton method to solve the nonlinear
                                        %     complex-valued system of equations directly
                                        % 12: Use self-programmed Newton method to solve the nonlinear
                                        %     complex-valued system of equations by splitting in real and imaginary part
                                        % 21: Use "fsolve" to solve the nonlinear complex-valued system
                                        %     of equations directly
                                        % 22: Use "fsolve" to solve the nonlinear complex-valued system
                                        %     of equations by splitting in real and imaginary part
                                        % if solution_method < 0:
                                        % 11 or 21 : Iterate the complex-valued solution vector
                                        % 12 or 22 : Iterate the real and complex parts of the
                                        %            solution vector separately
  % Used if solution_method > 0 & solver_method = 11 or 12:
  TolX = 1e-6;                          % if solver_method = 2: Error tolerance in iterate
  TolF = 1e-6;                          % if solver_method = 2: Error tolerance in function value
  itermax = 50;                         % if solver_method = 2: Maximum number of Newton iterations

  % Used if solution_method > 0 & solver_method = 21 or 22:
  TolX_nleq = 1e-6;                     % Error tolerance in iterate (default value: 1e-6)
  TolFun_nleq = 1e-6;                   % Error tolerance in function value (default value: 1e-6)
  MaxIter_nleq = 400;                   % Maximum number of Newton iterations (default value: 400)
  MaxFunEvals_nleq = [];                % Maximum number of function evaluations (default value: [])
  AutoScaling_nleq = 'off';             % default value: 'off'
  ComplexEqn_nleq = 'off';              % default value: 'off'
  Updating_nleq = 'off';                % Use Broyden updates to approximate the Jacobian (default value: 'off')
  TypicalX_nleq = [];                   % Typical order of the solution (default value : [])
  Jacobian_nleq = 'off';                % User-supplied Jacobian (default value: 'off')

  % Used if solution_method = 0
  urelax = 1;                           % Underrelaxation factor
  tol = 1e-6;                           % convergence tolerance
  max_iter = 50;                        % max  iterations

  % Set integration method
  integral_method = 0;                  % 0: Use integral, 1: Use trapezoidal rule
  RelTol_integral = 1e-6;               % if integral_method = 0: Relative error tolerance in integral solving (default value: 1e-6)
  AbsTol_integral = 1e-10;              % if integral_method = 0: Absolute error tolerance in integral solving (default value: 1e-10)
  nr = 300;                             % if integral_method = 1: number of discretization points in r-direction

  % Set computational case
  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)
    d = 38;                             % Interface depth (m) % Draft (m)
    tau = 0.1;                          % porosity ratio
    g = 9.81;                           % Gravity (m/s^2)
    gamma = 1.0;                        % discharge coefficient
    l = 0;                              % Azimuthal order
    X_kRE_min = 0.0025;
    X_kRE_max = 0.16;                   %%%%%% this is now k_0 ;
    nXsteps = 100;
    stepX = (X_kRE_max-X_kRE_min)/...
             nXsteps;
    X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
    stepX_min = stepX/2^6;
  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
    g = 9.81;                           % Gravity (m/s^2)
    gamma = 1.0;                        % discharge coefficient
    l = 0;                              % Azimuthal order
    X_kRE_min = 2.0653;%0.05
    X_kRE_max = 2.5914;%4.5             %%%%%% this is now k_0 ;
    nXsteps = 160;%40
    stepX = (X_kRE_max-X_kRE_min)/...
             nXsteps;
    X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
    nXsteps = 5;
    X_kRE=X_kRE(1:6);
    stepX_min = stepX/2^15;
  end

  % Initialize solution vectors
  if solver_method == 11 || solver_method == 21
    if abs(solution_method) <= 2
      solini = zeros(Q,1);
    elseif abs(solution_method) == 3
      solini = zeros(N+Q,1);
    elseif abs(solution_method) == 4
      solini = zeros(2*N+2*M+Q,1);
    end
  elseif solver_method == 12 || solver_method == 22
    if abs(solution_method) <= 2
      solini = zeros(2*Q,1);
    elseif abs(solution_method) == 3
      solini = zeros(2*(N+Q),1);
    elseif abs(solution_method) == 4
      solini = zeros(2*(2*N+2*M+Q),1);
    end
  end

  % Compute starting solution
  x = X_kRE_min;
  [converged,iter,resnorm,exitflag,sol,t,eta0,...
   A_vec,B_vec,C_vec,D_vec,E_vec,...
   AmA_vec,BmB_vec,FmC_vec,DmD_vec,WmE_vec] = ...
     get_solution(N,M,Q,RI,RE,h,d,tau,x,g,gamma,l,...
                  solution_method,...
                  solver_method,...
                  TolX,TolF,itermax,...
                  TolX_nleq,TolFun_nleq,MaxIter_nleq,MaxFunEvals_nleq,...
                  AutoScaling_nleq,ComplexEqn_nleq,Updating_nleq,...
                  TypicalX_nleq,Jacobian_nleq,...
                  urelax,tol,max_iter,...
                  integral_method,RelTol_integral,AbsTol_integral,nr,...
                  solini);

  if exitflag == 1 || converged == true
    start = true;
  else
    start = false;
  end

  % If start in X_kRE_min is successful, continue solution for complete X_kRE vector
  if start
    disp("Starting computations ...");
    disp(" ");

    T    = zeros(1,nXsteps+1);
    Eta0 = zeros(1,nXsteps+1);
    T(1) = t;
    Eta0(1) = eta0;

    Tfull(1) = t;
    Eta0full(1)= eta0;

    ns = 1;
    %dx = stepX;
    %Dx = stepX;

    for i = 2:nXsteps+1
      Xstart = X_kRE(i-1);
      Xend   = X_kRE(i);
      %x = Xstart + Dx;
      x = Xend;
      dx = stepX;
      sol0 = sol;
      finished = false;
      nsucc = 0;
      while ~finished && dx >= stepX_min
        [converged,iter,resnorm,exitflag,sol,t,eta0,...
         A_vec,B_vec,C_vec,D_vec,E_vec,...
         AmA_vec,BmB_vec,FmC_vec,DmD_vec,WmE_vec] = ... ...
           get_solution(N,M,Q,RI,RE,h,d,tau,x,g,gamma,l,...
                        solution_method,...
                        solver_method,...
                        TolX,TolF,itermax,...
                        TolX_nleq,TolFun_nleq,MaxIter_nleq,MaxFunEvals_nleq,...
                        AutoScaling_nleq,ComplexEqn_nleq,Updating_nleq,...
                        TypicalX_nleq,Jacobian_nleq,...
                        urelax,tol,max_iter,...
                        integral_method,RelTol_integral,AbsTol_integral,nr,...
                        sol0);
        if (solution_method <= 0 && ~converged)  || ...
           (solution_method > 0  && exitflag ~=1)
          x
          dx
          eta0
          if solution_method <= 0
            converged
            iter
          else
            exitflag
          end
          resnorm
          disp(" ");
        end
        if exitflag == 1 || converged == true
          ns = ns + 1;
          Tfull(ns) = t;
          Eta0full(ns) = eta0;
          if abs(Xend - x) < 1e-8
            finished = true;
            break
          end
          %Dx = min(stepX,2*dx);
          %Dx = dx;
          nsucc = nsucc + 1;
          if nsucc >= 3
            dx = min(Xend-x,2*dx);
            nsucc = 0;
          else
            dx = min(Xend-x,dx);
          end
          x = x + dx;
          sol0 = sol;
        else
          nsucc = 0;
          x = x - dx;
          dx = dx/2;
          %Dx = dx;
          x = x + dx;
          sol0 = sol0;
        end
      end
      x
      dx
      eta0
      disp(" ");

      T(i) = t;
      Eta0(i) = eta0;
    end

    % Write results to file ...
    AA = [T.',Eta0.'];
    % ... for prescribed grid points
    fid = fopen('data_T', 'w+');
    for i = 1:size(AA, 1)
      fprintf(fid, '%f ', AA(i,:));
      fprintf(fid, '\n');
    end
    fclose(fid);

    AA = [Tfull.',Eta0full.'];
    % ... for complete set of used grid points
    fid = fopen('data_T_full', 'w+');
    for i = 1:size(AA, 1)
      fprintf(fid, '%f ', AA(i,:));
      fprintf(fid, '\n');
    end
    fclose(fid);

    hold on
    plot(T,Eta0)
    plot(Tfull,Eta0full)
    hold off
    grid on
    size(Eta0)
    size(Eta0full)
  else
    disp("Start was not successful.")
  end

  toc
end

function [converged,iter,resnorm,exitflag,sol,xeta0,eta0,...
          A_vec,B_vec,C_vec,D_vec,E_vec,...
          AmA_vec,BmB_vec,FmC_vec,DmD_vec,WmE_vec] = ...
            get_solution(N,M,Q,RI,RE,h,d,tau,X_kRE,g,gamma,l,...
                         solution_method,...
                         solver_method,...
                         TolX,TolF,itermax,...
                         TolX_nleq,TolFun_nleq,MaxIter_nleq,MaxFunEvals_nleq,...
                         AutoScaling_nleq,ComplexEqn_nleq,Updating_nleq,...
                         TypicalX_nleq,Jacobian_nleq,...
                         urelax,tol,max_iter,...
                         integral_method,RelTol_integral,AbsTol_integral,nr,...
                         sol0);

converged = false;
iter = 0;
resnorm = 0;
exitflag = 0;

%clear all;
%close all;
%clc;
%tic;
% diary('iter_log.txt');
% diary on
% fprintf('Run started: %s\n', datestr(now));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%          Parameters from your setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%if solution_method == 1 || solution_method == 2 || ...
%   solution_method == 3 || solution_method == 4
%  first = true;
%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
%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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A_vec = zeros(N,length(X_kRE));
B_vec = zeros(M,length(X_kRE));
C_vec = zeros(M,length(X_kRE));
D_vec = zeros(N,length(X_kRE));
E_vec = zeros(Q,length(X_kRE));
AmA_vec = zeros(M,length(X_kRE));
BmB_vec = zeros(N,length(X_kRE));
FmC_vec = zeros(N,length(X_kRE));
DmD_vec = zeros(M,length(X_kRE));
WmE_vec = zeros(N,length(X_kRE));
 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


% % 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 ];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if solution_method < 0             % Fixed-point iterations

  Y = sol0;
  f = @(Y) wrapper(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
                   solution_method,solver_method,...
                   integral_method,RelTol_integral,AbsTol_integral,nr);
  for iter = 1:max_iter
    [Y_new,a_vec,b_vec,c_vec,d_vec,e_vec] = f(Y);
    diff_Y = max(abs(Y_new - Y));
    if diff_Y < tol
      converged = true;
      break
    end
    Y = Y_new;
  end
  sol = Y_new;
  %
  %Jac = JacMatrix(f,sol);
  %rho = norm(Jac)
  resnorm = diff_Y;

elseif solution_method == 0     % Ali

  %converged = false;
  %psi = zeros(Q,1);
  psi = sol0;
  %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;
      if integral_method == 0
        psi(q) = quadgk(integrand, 0, RI,  ...
                        'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
%        psi(q) = integral(integrand, 0, RI,  ...
%                          'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
      elseif integral_method == 1
        R = linspace(0,RI,nr);
        psi(q) = trapz(R,integrand(R));
      end
    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
      T_old = urelax*T + (1-urelax)*T_old;
    end
  end  % <<< end of iter loop
  sol = psi;
  resnorm = diff_T;

elseif solution_method > 0       % Newton-based formulations

  if solver_method == 11 || solver_method == 12
    f = @(Y) wrapper(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
                    solution_method,solver_method,...
                    integral_method,RelTol_integral,AbsTol_integral,nr);
    df = @(Y) JacMatrix(@(Z)f(Z),Y);
    Yold = sol0;
    iter = 0;
    while iter < itermax
      iter = iter + 1;
      [dfY,fY,a_vec,b_vec,c_vec,d_vec,e_vec] = df(Yold);
      dY = -dfY\fY;
      Ynew = Yold + dY;
      ErrorX = max(abs(dY));
      ErrorF = max(abs(fY));
      if ErrorX < TolX && ErrorF < TolF
        exitflag = 1;
        break
      end
      Yold = Ynew;
    end
    iter
    sol = Ynew;
    resnorm = ErrorF;
  elseif solver_method == 21 || solver_method == 22
    f = @(Y) wrapper(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
                     solution_method,solver_method,...
                     integral_method,RelTol_integral,AbsTol_integral,nr);
    df = [];
    if strcmp(Jacobian_nleq,'on')
      df = @(Y)JacMatrix(@(Z)f(Z),Y);
%      df = @(Y)jacobianest(@(Z)f(Z),Y);
%      df = @(Y)numgradient('wrapper',{Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
%                                      solution_method,...
%                                      integral_method,RelTol_integral,AbsTol_integral,nr});
    end
    options = optimset('Display','none','TolX',TolX_nleq,'TolFun',TolFun_nleq,...
                       'MaxIter',MaxIter_nleq,'MaxFunEvals',MaxFunEvals_nleq,...
                       'AutoScaling',AutoScaling_nleq,'ComplexEqn',ComplexEqn_nleq,...
                       'Updating',Updating_nleq,'TypicalX',TypicalX_nleq,...
                       'Jacobian',Jacobian_nleq);
    [sol,resnorm,exitflag,~,Jac_fsolve] = fsolve({f,df},sol0,options);
    [~,a_vec,b_vec,c_vec,d_vec,e_vec] = f(sol);
    resnorm = norm(resnorm);
    %Jac_fsolve
    %Jac = JacMatrix(@(Z)f(Z),sol);
    %Jac = jacobianest(@(Z)f(Z),sol)
    %Jac = numgradient('wrapper',{sol,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
    %                             solution_method,solver_method,...
    %                             integral_method,RelTol_integral,AbsTol_integral,nr})
  end
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);
 xeta0(idx) = omega(idx)^2*RE/g;

 A_vec(:,idx) = a_vec;
 B_vec(:,idx) = b_vec;
 C_vec(:,idx) = c_vec;
 D_vec(:,idx) = d_vec;
 E_vec(:,idx) = e_vec;

 AmA_vec(:,idx) = A*a_vec;
 BmB_vec(:,idx) = B*b_vec;
 FmC_vec(:,idx) = F*c_vec;
 DmD_vec(:,idx) = D*d_vec;
 WmE_vec(:,idx) = W*e_vec;

%if solution_method == 0
%  idx
%  converged
%  iter
%  if ~converged
%    eta0(idx) = NaN;
%  endif
%elseif solution_method == 1 || solution_method == 2 || ...
%       solution_method == 3 || solution_method == 4
%  idx
%  if exitflag == 1
%    exitflag,exitflag_check
%    norm(resnorm),norm(resnorm_check)
%  else
%    exitflag
%    norm(resnorm)
%  endif
%  if exitflag ~= 1
%    eta0(idx) = NaN;
%  endif
end
%eta0(idx)

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']);

%end

function [fY,a_vec,b_vec,c_vec,d_vec,e_vec] = ...
              wrapper(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
                      solution_method,solver_method,...
                      integral_method,RelTol_integral,AbsTol_integral,nr)

  if solver_method == 11 || solver_method == 21
    [fY,a_vec,b_vec,c_vec,d_vec,e_vec] = ...
              compute_fY(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
                         solution_method,...
                         integral_method,RelTol_integral,AbsTol_integral,nr);
  elseif solver_method == 12 || solver_method == 22
    Yc = Y(1:end/2) + 1i*Y(end/2+1:end);
    [fYc,a_vec,b_vec,c_vec,d_vec,e_vec] = ...
              compute_fY(Yc,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
                         solution_method,...
                         integral_method,RelTol_integral,AbsTol_integral,nr);
    fY = [real(fYc);imag(fYc)];
  end

end


function [fY,a_vec,b_vec,c_vec,d_vec,e_vec] = ...
              compute_fY(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
                         solution_method,...
                         integral_method,RelTol_integral,AbsTol_integral,nr)

  if solution_method == -1 || solution_method == 1
    psi = Y;
    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;
    a_vec = T(M+1:M+N);
    b_vec = T(1:M);
    c_vec = T(M+N+1:2*M+N);
    d_vec = T(2*M+N+1:2*M+2*N);
    e_vec = T(2*M+2*N+1:2*M+2*N+Q);
    fpsi = zeros(size(psi));
    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;
      if integral_method == 0
        fpsi(q) = quadgk(integrand, 0, RI,  ...
                         'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
%        fpsi(q) = integral(integrand, 0, RI,  ...
%                           'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
      elseif integral_method == 1
        R = linspace(0,RI,nr);
        fpsi(q) = trapz(R,integrand(R));
      end
    end
    if solution_method == -1
      fpsi = fpsi;
    elseif solution_method == 1
      fpsi = psi - fpsi;
    end
    fY = fpsi;
  elseif solution_method == -2 || solution_method == 2
    T = Y;
    Ur = U(1:2*M+2*N) - S(1:2*M+2*N,2*M+2*N+1:2*M+2*N+Q)*T;
    Tr = S(1:2*M+2*N,1:2*M+2*N)\Ur;
    a_vec = Tr(M+1:M+N);
    b_vec = Tr(1:M);
    c_vec = Tr(M+N+1:2*M+N);
    d_vec = Tr(2*M+N+1:2*M+2*N);
    e_vec = T;
    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;
      if integral_method == 0
        psi = quadgk(integrand, 0, RI,  ...
                     'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
%        psi = integral(integrand, 0, RI,  ...
%                       'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
      elseif integral_method == 1
        R = linspace(0,RI,nr);
        psi = trapz(R,integrand(R));
      end
      U(2*M+2*N+q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi;
    end
    if solution_method == -2
      Unr = U(2*M+2*N+1:2*M+2*N+Q) - S(2*M+2*N+1:2*M+2*N+Q,1:2*M+2*N)*Tr;
      Tnr = S(2*M+2*N+1:2*M+2*N+Q,2*M+2*N+1:2*M+2*N+Q)\Unr;
      fT = Tnr;
    elseif solution_method == 2
      fT = S(2*M+2*N+1:2*M+2*N+Q,1:2*M+2*N)*Tr(1:2*M+2*N) + ...
           S(2*M+2*N+1:2*M+2*N+Q,2*M+2*N+1:2*M+2*N+Q)*T - ...
           U(2*M+2*N+1:2*M+2*N+Q);
    end
    fY = fT;
  elseif solution_method == -3 || solution_method == 3
    T = Y;
    Ur = U(1:2*M+N) - S(1:2*M+N,2*M+N+1:2*M+2*N+Q)*T;
    Tr = S(1:2*M+N,1:2*M+N)\Ur;
    a_vec = Tr(M+1:M+N);
    b_vec = Tr(1:M);
    c_vec = Tr(M+N+1:2*M+N);
    d_vec = T(1:N);
    e_vec = T(N+1:N+Q);
    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;
      if integral_method == 0
        psi = quadgk(integrand, 0, RI,  ...
                     'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
%        psi = integral(integrand, 0, RI,  ...
%                      'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
      elseif integral_method == 1
        R = linspace(0,RI,nr);
        psi = trapz(R,integrand(R));
      end
      U(2*M+2*N+q) = -(1i*b1/omega(idx)) * (2/(RI^2 * dbesselj(l, chi(q)*RI))) * psi;
    end
    if solution_method == -3
      Unr = U(2*M+N+1:2*M+2*N+Q) - S(2*M+N+1:2*M+2*N+Q,1:2*M+N)*Tr;
      Tnr = S(2*M+N+1:2*M+2*N+Q,2*M+N+1:2*M+2*N+Q)\Unr;
      fT = Tnr;
    elseif solution_method == 3
      fT = S(2*M+N+1:2*M+2*N+Q,1:2*M+N)*Tr + ...
           S(2*M+N+1:2*M+2*N+Q,2*M+N+1:2*M+2*N+Q)*T - ...
           U(2*M+N+1:2*M+2*N+Q);
    end
    fY = fT;
  elseif solution_method == -4 || solution_method == 4
    T = Y;
    a_vec = T(M+1:M+N);
    b_vec = T(1:M);
    c_vec = T(M+N+1 : 2*M+N);
    d_vec = T(2*M+N+1:2*M+2*N);
    e_vec = T(2*M+2*N+1:end);
    psi = zeros(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;
      if integral_method == 0
        psi(q) = quadgk(integrand, 0, RI,  ...
                        'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
%        psi(q) = integral(integrand, 0, RI,  ...
%                          'RelTol',RelTol_integral,'AbsTol',AbsTol_integral);
      elseif integral_method == 1
        R = linspace(0,RI,nr);
        psi(q) = trapz(R,integrand(R));
      end
    end
    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
    if solution_method == -4
      fT = S\U;
    elseif solution_method == 4
      fT = S*T - U;
    end
    fY = fT;
  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

function [Jac,fY,a_vec,b_vec,c_vec,d_vec,e_vec] = JacMatrix(f,Y)

  % Forward differencing

  [fY,a_vec,b_vec,c_vec,d_vec,e_vec] = f(Y);

  dimY = numel(Y);
  Jac = zeros(dimY);
  for j = 1:dimY
    Ysave = Y(j);
    % Choose deltaY as to keep increment in the direction of Y(j)
    if Ysave == 0
      deltaY = sqrt(eps*1e-5);
    else
      deltaY = sqrt(eps)*sqrt(max(1e-5,abs(Ysave)))*Ysave/abs(Ysave);
    end
    %% Stepsize control by Hairer/Wanner in ODE integrator RADAU5
    %deltaY = sqrt(eps*max(1e-5,abs(Ysave)));
    %% Stepsize control in solver hybrd for systems of nonlinear equations
    %deltaY = sqrt(eps)*max(abs(Ysave),1)
    Y(j) = Y(j) + deltaY;
    fYvar = f(Y);
    Y(j) = Ysave;
    Jac(:,j) = (fYvar - fY)/deltaY;
  end

  % Central differencing
  %
  %dimY = numel(Y);
  %Jac = zeros(dimY);
  %for j = 1:dimY
  %  Ysave = Y(j);
  %  if Ysave == 0
  %    deltaY = sqrt(eps*1e-5);
  %  else
  %    deltaY = sqrt(eps)*sqrt(max(1e-5,abs(Ysave)))*Ysave/abs(Ysave);
  %  endif
  %  %deltaY = sqrt(eps*max(1e-5,abs(Ysave)));
  %  %deltaY = sqrt(eps)*max(abs(Ysave),1)
  %  Y(j) = Y(j) + deltaY;
  %  fYvarp = compute_fY(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
  %                      solution_method,...
  %                      integral_method,RelTol_integral,AbsTol_integral,nr);
  %  Y(j) = Ysave;
  %  Y(j) = Y(j) - deltaY;
  %  fYvarm = compute_fY(Y,S,U,M,N,Q,Z0d,wk,RI,l,Znd,k,chi,b1,omega,idx,...
  %                      solution_method,...
  %                      integral_method,RelTol_integral,AbsTol_integral,nr);
  %  Y(j) = Ysave;
  %  for i = 1:dimY
  %    Jac(i,j) = (fYvarp(i) - fYvarm(i))/(2*deltaY);
  %  end
  %end

end

EDIT:

At the beginning of the code, change

  % Initialize solution vectors
  if solver_method == 11 || solver_method == 21
    if abs(solution_method) <= 2
      solini = zeros(Q,1);
    elseif abs(solution_method) == 3
      solini = zeros(N+Q,1);
    elseif abs(solution_method) == 4
      solini = zeros(2*N+2*M+Q,1);
    end
  elseif solver_method == 12 || solver_method == 22
    if abs(solution_method) <= 2
      solini = zeros(2*Q,1);
    elseif abs(solution_method) == 3
      solini = zeros(2*(N+Q),1);
    elseif abs(solution_method) == 4
      solini = zeros(2*(2*N+2*M+Q),1);
    end
  end

to

  % Initialize solution vectors
  if solver_method == 11 || solver_method == 21
    if solution_method == 0
      solini = zeros(Q,1);
    elseif abs(solution_method) == 1 || abs(solution_method) == 2
      solini = zeros(Q,1);
    elseif abs(solution_method) == 3
      solini = zeros(N+Q,1);
    elseif abs(solution_method) == 4
      solini = zeros(2*N+2*M+Q,1);
    end
  elseif solver_method == 12 || solver_method == 22
    if solution_method == 0
      solini = zeros(Q,1);
    elseif abs(solution_method) == 1 || abs(solution_method) == 2
      solini = zeros(2*Q,1);
    elseif abs(solution_method) == 3
      solini = zeros(2*(N+Q),1);
    elseif abs(solution_method) == 4
      solini = zeros(2*(2*N+2*M+Q),1);
    end
  end