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

[Javeria]

@Torsten thanks for your guidance . I checked the code as per your instructions and yes it takes few seconds but as i increased the N=M=Q=20,50,80 so one because i also wants to check that how many terms are enough for my series convergence i.e how many terms is sufficient and after that number if i increased this number there is no effect on the solution. but by doing so the code takes time. I just want to see the results if we go with the traphizoidal rule for the integral i.e

nr = 300; % if integral_method = 1: number of discretization points in r-direction

"However the code you provided have so much solvers" so i just want this one solver for my problem to test it by varying the N=M=Q. could you please help me by just dropping the script which only include this one solver

thanks

[Torsten]

First change

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

to

%    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);
    options = optimset('Display','none','TolX',TolX_nleq,'TolFun',TolFun_nleq,...
                       'MaxIter',MaxIter_nleq,'MaxFunEvals',MaxFunEvals_nleq); 

Now you should also be able to work with MATLAB's "fsolve".

Test this with a small runcase for the following settings:

  solution_method = 2;
  solver_method = 22;
  % 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
  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 = 0.05
  X_kRE_max = 4.5                     %%%%%% this is now k_0 ;
  nXsteps = 40
  stepX = (X_kRE_max-X_kRE_min)/...
           nXsteps;
  X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
  stepX_min = stepX/2^15;  

Since it's not certain whether there is one solver which works for all your applications, I don't want to shorten the code. But in essence, there are only 4 settings that you can play with:

solution_method (recommended is = 1 or = 2)

solver_method (recommended is = 12 or - if "fsolve" works with the example above - = 22)

integral_method (= 0 if you want to use the adaptive MATLAB integrator, = 1 if you want to use the trapezoidal rule). If you choose integral_method = 1, you additionally have to set "nr" to the number of evaluation points in the interval [0 RI] for the trapezoidal method.

If it's too confusing for you, you can write a script and call "main" with your settings:

  solution_method = ...;
  solver_method = ...;
  integral_method = ...
  nr = ...;
  M = ...;
  N = ...;
  Q = ...;
  main(solution_method, solver_method, integral_method, nr, M, N, Q)

Of course, you have to comment out the respective settings in "main" before for that your values are used and not the ones defined in "main".

You should write intermediate results for eta0 to file, not only at the end when eta0 is computed for all values of XkRE (as it is done in the actual code). Even saving all relevant parameters (T,psi,...) for a restart at the value for X_kRE when the solver potentially crahed is advisable.

But be prepared that running cases for big values for M/N/Q (especially Q) will take a long time.

Shouldn't testing

 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; 

suffice for this purpose ?

[Torsten]

@Javeria

Since I don't use MATLAB at the moment, I'd be interested in runtime information for the cases you solve to compare them with octave performance.

Is it possible for you to provide this information together with the values for solution_method, solver_method, integration_method, nr and the case parameters (N,M,Q,RI,RE,h,d,tau,gamma,l) you used ?

[Javeria]

@Torsten i just give M=N=Q=50 and tested for the small case

nXsteps = 5;

X_kRE=X_kRE(1:6);

and i have the following .

Elapsed time is 8316.962963 seconds.

[Javeria]

@Torsten i run your code with the following setup

main()
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 = 40;                             % Number of N modes
    M = 40;                             % Number of M modes
    Q = 40;                             % 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 =0.05; %2.0653;%0.05
    X_kRE_max = 4.5;%2.5914;%4.5             %%%%%% this is now k_0 ;
     % nXsteps = 5;
       nXsteps = 40;%160;%40
     stepX = (X_kRE_max-X_kRE_min)/...
             nXsteps;
     X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
    
     %X_kRE=X_kRE(1:6);
    stepX_min = stepX/2^15;
  end

and i get these results attached as png and the time for this simulation is

Elapsed time is 23803.745739 seconds.

>>

[Torsten]

Thank you for the infos.

I also ran the code under octave and got this result for the M = N = Q = 20 case (with no other settings changed except solver_method = 22 and Updating_nleq = 'on' ) within appr. 6 hours. Unfortunately, the "Updating" option is not available under MATLAB - it speeds things up quite significantly, but makes the computation slightly more unstable.

Here is a comparison with the M = N = Q = 10 case:

A1 = readmatrix('data_T_20_20_20_fein_full.txt');
A2 = readmatrix('data_T_10_10_10_fein.txt');
hold on
plot(A1(:,1),A1(:,2))
plot(A2(:,1),A2(:,2))
hold off
grid on

●

I'm curious and take your case to see how long it takes under octave.

[Javeria]

@Torsten yes so which solver i need to use to solve the problem and take less time.

[Torsten]

At present, for larger problems, I use

solution_method = 2
solver_method = 22
Update_nleq = 'on'
integration_method = 0

with octave.

MATLAB does not have the "Update" option which makes the code significantly faster (but a bit more unstable). And I cannot test MATLAB's "fsolve" and see whether you won't come into trouble with its "exitflag". In the part of the code I wrote, exitflag = 1 indicates that the computation can progress to the next grid point in X_kRE. If exitflag ~= 1, the code tries to insert intermediate points in X_kRE to see whether it can now successfully continue. This can cost quite a lot of time for larger problems and can even make the code stop. And I don't know whether MATLAB's "fsolve" reliably returns exitflag = 1 as the "fsolve" solver under "octave" does (the two solvers are not identical although both share the same name).

Further, until now, I didn't test whether integral_method = 1 makes the code faster. This would also be an option you could test.

So as already adviced above:

Try whether "fsolve" without the Update option (solver_method = 22 instead of solver_method = 12) is still faster than my primitive Newton method code and doesn't make the stepsize being reduced because it doesn't return exitflag = 1. And test integral_method = 1 with an appropriate value for "nr" instead of integral_method = 0.

Or at least test octave instead of MATLAB.

Writing a numerical code takes some time. But most of the time is spent with testing and finding the best solution method.

[Torsten]

@Javeria

Does your computer and your MATLAB licence admit parallel computing ?

Until now, we use the solution for eta0 in point X_kRE(i-1) as initial guess for the solution in point X_kRE(i) because it stabilizes the computation. So this version of the code would not be parallizable.

But if the solver succeeded also if you could solve for eta(i) independent of what eta0 was in X_kRE(i-1), you could compute in parallel.

[Javeria]

@Torstenlicense('test','Distrib_Computing_Toolbox')

ans =

1

i tested my MATLAB

[Torsten]

First test whether the code works if you don't use the solution of the last step for a restart without parallel computing:

Replace

% 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
  

by

  fid = fopen('data.txt', 'w+');
  T    = zeros(1,nXsteps+1);
  Eta0 = zeros(1,nXsteps+1);
  for i = 1:nXsteps+1
    x = X_kRE(i);
   [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); 
    T(i) = t;
    Eta0(i) = eta0;
    fprintf(fid, '%f %f \n ',t,eta0);
  end
  fclose(fid);
  plot(T,Eta0)
  grid on
  

and see if it works for the case

  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 = 0.05
  X_kRE_max = 4.5                     %%%%%% this is now k_0 ;
  nXsteps = 40
  stepX = (X_kRE_max-X_kRE_min)/...
           nXsteps;
  X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);
  

If yes, try parallelizing the loop with "parfor". I have no experience with parallel computing - so if you have problems, you will need to open a new question.

Answer 2

[William Rose]

@Javeria, @Torsten,

I am learning from your discussion. I also learned from the comments of @Sam Chak and @David Goodmanson. I like Sam's reminder to focus on the issues that most pertain to the original research question.

@Javeria, I strongly recommend that you accept @Torsten's answer.

[Javeria]

@Torsten@David Goodmanson@Sam Chak@Walter Roberson Thanks to all of you for your suggestions and guidance. However using the above code it also increase computational time for me as i am also testing the convergence for the N modes. I did the following correction because i note that the integral blows up for the r values. and then that's why i have NaN for some r. Indeed in my analytical formulation the domain which we have is from 0 to RI so that's why i take then integration of that domain from 0 to RI but in numerical formulation this cause instabilities for small RI. so to prevent from it i did the following changes and it give me the solution but as when i increased the number of eigen mode i.e N =80 upt0 60 it works then again my code give me some numerical instabilities. Any comments and any suggestions about this instabilities i want a robust code which works for all N and for all k0. below is my 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 = 20;                    % Number of N modes
M = 20;                    % Number of M modes
Q = 20;                    % Number of Q modes
% RI = 34.5;                 % Inner radius (m)
% ratio = 47.5/34.5;         % Ratio R_E / R_I
% % ratio = 2.0;
% RE = ratio*RI;             % outer radius (m)
% h = 200/34.5 * RI;         % Water depth (m)
% d = 38/34.5 * RI;              % Interface depth (m)
RI = 0.345;
RE =0.475;
h = 2.0;
d= 0.38;
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% T1 = linspace(5, 35, 50);           % Define wave period range [s]
% fprintf('Check T1 immediately after linspace: %d\n', length(T1));
% omega = 2 * pi ./ T1;               % Angular frequency [rad/s]
kh = linspace(0.1,20,50);
omega = zeros(size(kh));         % Preallocate omega
eta0 = zeros(size(kh));          % Preallocate eta0
% X_kRE = zeros(size(omega));        % Preallocate k_0 array
 % % X_kRE = linspace(0.0025, 0.16,100);         %%%%%% this is now k_0 
% kh_range = linspace(0.01, 10, 100);   % Covering shallow to deep regimes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% eta0 = zeros(size(X_kRE));  % Preallocate
% for ii = 1:length(omega)
%     dispersion_eq = @(k) omega(ii)^2 - g * k * tanh(k * h);
%     X_kRE(ii) = fzero(dispersion_eq, [1e-6, 10]);
% end
for idx = 1:length(kh)
    % if I_given_wavenumber == 1
        wk = kh(idx)/h;                % dimensional wavenumber  
        omega(idx) = sqrt(g * wk * tanh(wk * h)); % dispersion relation
       % fun_alpha = @(x) omega(idx)^2 + g*x.*tan(x*h);     %dispersion equation for evanescent modes

% % 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) debug_integrand(r, N, Q, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi, q);
r_min = max(1e-3 * RI, 1e-6);   % adaptive lower bound


        psi(q) = integral(integrand, 0,r_min, '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
% 
% % % ---------- Summary of iterations for T in [5, 35] s ----------
% % T_all = 2*pi./omega(:);
% % mask  = (T_all >= 5) & (T_all <= 35);
% % [Ts, order] = sort(T_all(mask), 'descend');   % sort by T (optional)
% % idx_list    = find(mask);
% % idx_list    = idx_list(order);
% % iters_list  = iters_used(mask);
% % iters_list  = iters_list(order);
% 
% % fprintf('\nIterations for T in [5, 35] s (sorted by T):\n');
% % fprintf('  idx   T [s]   iters\n');
% % for i = 1:numel(idx_list)
% %     fprintf(' %4d  %7.3f     %2d\n', idx_list(i), Ts(i), iters_list(i));
% % end
% 
% % % (optional) write the same table to a CSV file
% % tbl = table(idx_list, Ts, iters_list, eta0(idx_list)', ...
% %     'VariableNames', {'idx','T','iters','eta0'});
% % writetable(tbl, 'iters_5to35.csv');
% % fprintf('Wrote iters_5to35.csv with %d rows.\n', height(tbl));
% 
% % % Show which index is the resonance and how many iterations it used
% % [~, i_peak] = max(eta0);
% % fprintf('\nRES peak at idx %d (T=%6.3f s): iters=%d\n', ...
% %         i_peak, 2*pi/omega(i_peak), iters_used(i_peak));
% % 
% % % (optional) print a small window around resonance
% % win = 5;  % ±5 indices around the peak
% % i0  = max(1, i_peak - win);
% % i1  = min(length(X_kRE), i_peak + win);
% % fprintf('\nIterations near resonance (idx %d):\n', i_peak);
% % fprintf('  idx   T [s]   iters\n');
% % for j = i0:i1
% %     fprintf(' %4d  %7.3f     %2d\n', j, 2*pi/omega(j), iters_used(j));
% % end
% 
% % % (optional) quick plot: iterations vs T
% % figure(1); clf;
% % plot(T_all, iters_used, '-o'); grid on;
% % xlabel('T [s]'); ylabel('iterations'); title('Nonlinear iterations vs T');
% % xlim([5 35]);
% 
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % %                   Load CSV files
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % Load CSV files HERE (after the loop ends)
%     % % data_chen = load('dr_chen_no_dissp.csv');
% % % data_liu  = load('dr_liu.csv');
%      % data_exp  = load('experimental.csv');
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % %              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 
% fprintf('Length of T1: %d\n', length(T1));
% fprintf('Length of eta0: %d\n', length(eta0));
% plot(T1, eta0, 'k', 'LineWidth', 1.5); %%% t
plot(kh, eta0, 'k', 'LineWidth', 1.5);
%     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');

Warning: Ignoring extra legend entries.

% 
  grid on;
  xlim([0 20]);  % Match the reference plot
  ylim([0 4.5]);  % Optional: based on expected peak height

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

Time consuming =  25.9534 s

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi)

    % --- Allow vector r during integration ---
    if numel(r) > 1
        v_D_val = zeros(size(r));
        for i = 1:numel(r)
            v_D_val(i) = v_D(N, Q, r(i), Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi);
        end
        return
    end

    % === Term 1: n = 0 ===
    denom0 = dbesselj(l, wk * RI);
    if ~isfinite(denom0) || abs(denom0) < 1e-10
        fprintf('⚠️ Problem in term1: dbesselj = %.4e\n', denom0);
        keyboard
    end
    term1 = (d_vec(1) * Z0d * besselj(l, wk * r)) / denom0;

    % === sum1: Evanescent modes (n ≥ 2) ===
    sum1 = 0;
    for nidx = 2:N
        arg_n = k(nidx) * RI;
        denom_i = dbesseli(l, arg_n);

        if ~isfinite(denom_i) || abs(denom_i) < 1e-10
            fprintf('⚠️ Problem in sum1 at nidx = %d\n', nidx);
            fprintf('dbesseli(l, k*RI) = %.4e\n', denom_i);
            keyboard
        end

        sum1 = sum1 + d_vec(nidx) * Znd(nidx) * besseli(l, k(nidx) * r) / denom_i;
    end

    % === sum2: Radial modes (q = 1..Q) ===
    sum2 = 0;
    for qidx = 1:Q
        bj  = besselj(l, chi(qidx) * r);
        dbj = dbesselj(l, chi(qidx) * RI);

        if ~isfinite(bj) || ~isfinite(dbj) || abs(dbj) < 1e-10
            fprintf('\n⚠️ Issue in sum2 (radial)\n');
            fprintf('qidx = %d | r = %.4e | chi = %.4e\n', qidx, r, chi(qidx));
            fprintf('besselj = %.4e | dbesselj = %.4e\n', bj, dbj);
            keyboard
        end

        sum2 = sum2 + e_vec(qidx) * chi(qidx) * bj / dbj;
    end

    % === Total value ===
    v_D_val = term1 + sum1 + sum2;

    % === Final check ===
    if ~isfinite(v_D_val)
        fprintf('\n‼️ Final NaN/Inf in v_D_val\n');
        fprintf('r = %.4e | term1 = %.4e | sum1 = %.4e | sum2 = %.4e\n', ...
            r, term1, sum1, sum2);
        keyboard
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function val = debug_integrand(r, N, Q, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi, q)
    vD = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi);
    J = besselj(l, chi(q) * r);
    val = abs(vD) .* vD .* J .* r;

    if any(~isfinite(val))
        fprintf('\n‼️ NaN/Inf in integrand at r = %.4e\n', r);
        fprintf('v_D = %.4e, besselj = %.4e, r = %.4e\n', vD, J, r);
        keyboard
    end
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 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 = 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 = 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 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

[Torsten]

@Javeria

I don't understand why you think

       psi(q) = integral(integrand, 0,r_min, 'AbsTol',1e-8, 'RelTol',1e-6);   

is approximately

       psi(q) = integral(integrand, 0,RI, 'AbsTol',1e-8, 'RelTol',1e-6);   

?

Don't you want to replace

       psi(q) = integral(integrand, 0,RI, 'AbsTol',1e-8, 'RelTol',1e-6);

by

       psi(q) = integral(integrand, r_min ,RI, 'AbsTol',1e-8, 'RelTol',1e-6);  

here ?

Further, you can avoid "integral" calling your function "debug_integrand" with r being a vector with more than one element by using 'Arrayvalued',true as further option when calling "integral":

       value_integral = integral(fun,a,b,'AbsTol',...,'RelTol',...,'ArrayValued',true);

This would turn your recursive calls of v_D in this part

       % --- Allow vector r during integration ---
       if numel(r) > 1
           v_D_val = zeros(size(r));
           for i = 1:numel(r)
               v_D_val(i) = v_D(N, Q, r(i), Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi);
           end
           return
       end

unnecessary.

Did you validate your results with the results of one of the Newton solvers ?

Here is a plot of your method against the Newton solver in the critical region:

 A = load('Results_Newton.txt');
 B = load('Results_Iteration.txt');
 hold on
 plot(A(:,1),A(:,2))
 plot(B(:,1),B(:,2))
 hold off
 grid on

●

Did you check how large M, N and Q have to be for the parameter set you use at the moment such that your series can be considered as converged ?

[Javeria]

@Torsten Thanks. Before i tested your solvers but when i increased the N number of turncation of series then that solver take too much time. Yes i need to check the convergence of number of turncation also and we can tested it as that after further increasing this number N then we have no change on our solution so then that number is enough in my large parameters and fixed point method i take N=90 is enough.

[Torsten]

I used integral_method = 1, nr = 160 , M = N = Q = 60 and the MATLAB version of the solver "nleq1" given at the below link and it took about 11 hours to run under "octave" for 160 output points. You can multiply the 11 hours by approximately 2/3 for MATLAB, I guess, so that you arrive at approximately 7 hours. The solution for eta0 stabilized at these values.

https://www.zib.de/codelib/NewtonLib/

Also using "fsolve" gave acceptable run times of about 15 hours under "octave", thus approximately 10 hours with MATLAB.

[Javeria]

@Torsten CAN YOU PLEASE SHARE THE UPDATED CODE. THANKS

[Torsten]

I think the settings in the attached code are best-possible for your case.

You will have to download the nleq1 code from the website because I'm not allowed to distribute it. Only copy the nleq1.m file - the other files are example files that are unnecessary.

[Javeria]

@Torsten OKAY THANKS

[Torsten]

Are you angry or what do your upper case letters mean ?

[Javeria]

@Torsten no i am not angry your suggestions are quite helpful for me and i am trying to get my required solution my capslock was on that time that why upper case letters was typed.

[Torsten]

I made some numerical experiments to find the reason why your fixed-point iteration works for one parameter set while it doesn't for the other.

Let's start simple to see the idea.

Let's compare the recursion x_(n+1) = 2*x_n with the recursion x_(n+1) = 0.5*x_n and an initial value x_0 ~=0. Then the first recursion diverges while the second converges to 0. Why is this so ?

A recursion is given by an iteration function f in the sense that x_(n+1) = f(x_n). For the above two cases, f(x) = 2*x for the first recursion and f(x) = 0.5*x for the second. A sufficient (and "almost" necessary) condition for convergence is |f'(x)| < 1 near the fixed point. This will ensure that if you start with x0 near the fixed point, the fixed point will "attract" the iterates.

If you generalize this fact to higher dimensions, the eigenvalues of the Jacobian matrix of f in the fixed point are the suitable indicator whether iteration will converge or not. So I computed the solution T* for the values of X_kRE with Newton iteration and deduced from it the Jacobian of the iteration function f at T*.

It turned out that in case of the parameter set where you got convergence, for all values of X_kRE, the Jacobian matrix of f in the T*'s only had eigenvalues of absolute values < 1 when you use the following setup:

    M = 20;                             % Number of M modes
    N = 20;                             % Number of N 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);

In case of the parameter set where you encounter divergence, the Jacobian at T* has eigenvalues in absolute value > 1 starting from X_kRE(15) when using the below setup (I didn't check when it ended) - and the fixed-point iteration started to diverge from X_kRE(17):

    M = 20;                             % Number of M modes
    N = 20;                             % Number of N modes
    Q = 20;                             % Number of Q modes
    RI = 0.345;                           % Inner radius (m)
    RE = 0.475;                           % Outer radius (m)
    h = 2.0;                            % Water depth (m) % Interface depth (m)
    d = 0.38;                            % 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 = 0.1;%0.05;%2.0653;%0.05;
    X_kRE_max = 4.2105;%20;%4.5;%2.5914;%4.5;             %%%%%% this is now k_0 ;
    nXsteps = 120;%40;%160;%40;
    stepX = (X_kRE_max-X_kRE_min)/...
             nXsteps;
    X_kRE = linspace(X_kRE_min,X_kRE_max,nXsteps+1);    

I hope I could convince you that using the fixed-point iteration approach independent of the parameter set doesn't make sense. It would be necessary to change the iteration function f in T_(n+1) = f(T_n) to have only eigenvalues of absolute value < 1 near the fixed points T*, but I have no idea how this could be achieved.

Answer 3

[idris]

Truly slow @Torsten, attached is the result

[Torsten]

We are still in search of a reliable and fast solution for the case where fixed point iteration cannot work. Do you have a suggestion ?

If it had been obvious right from the beginning that the project would become computationally intensive, I had chosen C or Fortran. Most probably only 1/10 of the computational time needed compared to MATLAB.

[Javeria]

@Torsten To solve nonlinear equations, there is not a general rule but each pb could apply different methods. In your case of quadratic damping, have you tried with an estimation of damping (linear or nonlinear) as the input? (these comments is from my professor)

[Torsten]

The software codes provided to solve the nonlinear part of your system of equations (my own simple Newton method, "fsolve" from MATLAB and "nleq1" from ZIB) are general-purpose codes. That means that they are not adopted in any way to solve your special problem for the Q modes. And for this kind of codes, there is nothing else but the formulation of the system of equations itself (thus the f in f(x) = 0) that you can use to influence their performance.

I don't know exactly what you mean by "estimation of damping" as input. I just took your equations and handed them to the solver without backgound knowledge about the physics of your problem. But if you intend to use "estimation of damping" as a means to influence the treatment of the nonlinearity and/or to accelerate convergence, my guess is that it won't be possible in the framework of these general purpose codes. To take advantage of special structures or features of your system to be solved, I guess you will have to implement it in your own nonlinear solver.

[idris]

@Javeria, is the problem solved now, as am seeeing you focusing on @Torsten's?

[Javeria]

@idris no, i did not solve my problem

[Javeria]

@Torsten see the below pdf and the equation number 24. They have the same non linear integral offcourse their expressions are changed but the non linear condition is same. They mentioned the method that how to solve it.

Answer 4

[Javeria]

@Torsten i have the following code for the numerical integration but i am not sure whether it will work for my case and then how i can use it for my case.

clear

% Consider f(x) = cos(x) * exp(sin(x)), i.e. d/dx [exp(sin(x))] 
f = @(x) cos(x) .* exp(sin(x));
a = 0.6; b = 2.2;


% The analytical evaluation is exp(sin(b)) - exp(sin(a))
int_analytical = exp(sin(b)) - exp(sin(a));

% Set the number of grid points
N = 101;


% Choose a method 'T', 'S', or 'G'. Need N odd for 'S', but N can be odd or
% even for 'T' and 'G'
method = 'T';

[x,w] = get_quadrature_weights(a, b, N, method);

int_numerical = sum(w.*f(x));

disp(['Analytical evaluation: ' num2str(int_analytical, 10)])

Analytical evaluation: 0.4857117348

disp(['Numerical evaluation: ' num2str(int_numerical, 10)])

Numerical evaluation: 0.4856852319

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,w] = get_quadrature_weights(a, b, N, method)
% get quadrature nodes (x) and weights (w) with N points over the 
% interval a < x < b. The method can be 'T' (trapezium rule), 'S' 
% (Simpson's rule) or 'G' (Gauss-Legendre quadrature).

% Trapezium rule: the spacing between points is equal
% Simpson's rule: N needs to be odd

if strcmp(method, 'T')
    dx = (b-a)/(N-1);
    x = linspace(a,b,N);
    w = dx/2 * [1, 2*ones(1,N-2), 1];
    
elseif strcmp(method, 'S')
    if mod(N,2) == 0
        error('N must be odd to use Simpson''s rule')
        return
    end
    
    dx = (b-a)/(N-1);
    x = linspace(a,b,N);
    n = (N-1)/2;
    A = [4*ones(1,n-1); 2*ones(1,n-1)];
    w = dx/3 * [1, A(:).', 4, 1];
    
elseif strcmp(method, 'G')
    [x,w] = lgwt(N,a,b);
else 
    error('Method needs to be "T", "S" or "G"')
    x = NaN; w = NaN;
end

end

% ------------------------------------------------------------------

function [x,w]=lgwt(N,a,b)
% lgwt.m
%
% This script is for computing definite integrals using Legendre-Gauss 
% Quadrature. Computes the Legendre-Gauss nodes and weights  on an interval
% [a,b] with truncation order N
%
% Suppose you have a continuous function f(x) which is defined on [a,b]
% which you can evaluate at any x in [a,b]. Simply evaluate it at all of
% the values contained in the x vector to obtain a vector f. Then compute
% the definite integral using sum(f.*w);
%
% Written by Greg von Winckel - 02/25/2004
N=N-1;
N1=N+1; N2=N+2;
xu=linspace(-1,1,N1)';
% Initial guess
y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);
% Legendre-Gauss Vandermonde Matrix
L=zeros(N1,N2);
% Derivative of LGVM
Lp=zeros(N1,N2);
% Compute the zeros of the N+1 Legendre Polynomial
% using the recursion relation and the Newton-Raphson method
y0=2;
% Iterate until new points are uniformly within epsilon of old points
while max(abs(y-y0))>eps
    
    
    L(:,1)=1;
    Lp(:,1)=0;
    
    L(:,2)=y;
    Lp(:,2)=1;
    
    for k=2:N1
        L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
    end
 
    Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);   
    
    y0=y;
    y=y0-L(:,N2)./Lp;
    
end
% Linear map from[-1,1] to [a,b]
x=(a*(1-y)+b*(1+y))/2;      
% Compute the weights
w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
end

[Javeria]

@Torsten i also used the following code where i have the scale mode (1:100) i.e for small parameters setup and i have the NaN for eta0

clear all;
close all;
clc;
tic;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%          Parameters - 1:100 SCALE MODEL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 10;                    % Number of N modes
M = 10;                    % Number of M modes
Q = 10;                    % Number of Q modes

% === SCALED GEOMETRIC PARAMETERS (divided by 100) ===
RI = 0.345;                % Inner radius (m) [34.5/100]
ratio = 47.5/34.5;         % Ratio R_E / R_I (dimensionless - unchanged)
RE = ratio*RI;             % outer radius (m) [0.4750 m]
h = 200/34.5 * RI;         % Water depth (m) [2.0 m]
d = 38/34.5 * RI;          % Interface depth (m) [0.38 m]

b = h-d;                   % Distance from cylinder bottom to seabed (m) [1.62 m]

% === PHYSICAL CONSTANTS (unchanged for Froude scaling) ===
g = 9.81;                  % Gravity (m/s^2) - UNCHANGED
tau = 0.2;                 % porosity ratio - UNCHANGED (dimensionless)
gamma = 1.0;               % discharge coefficient - UNCHANGED (dimensionless)

b1    = (1 - tau) / (2 * gamma * tau^2);   % nonlinear coeff - UNCHANGED
tol      = 1e-4;           % convergence tolerance
max_iter = 100;            % max iterations
l = 0;                     % Azimuthal order

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === MOLIN'S DISCRETIZATION SETUP (DO THIS ONCE) ===
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Npt = 200;                           % Number of radial segments (Molin used 200)
R_seg = linspace(0, RI, Npt+1);      % Segment boundaries
R_mid = (R_seg(1:end-1) + R_seg(2:end)) / 2;  % Midpoints for evaluation

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === SCALED TIME PERIOD RANGE (divided by √100 = 10) ===
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T1 = linspace(0.5, 3.5, 10);        % Scaled period range [s] [5/10 to 35/10]
fprintf('Check T1 immediately after linspace: %d\n', length(T1));

Check T1 immediately after linspace: 10

omega = 2 * pi ./ T1;               % Angular frequency [rad/s]
X_kRE = zeros(size(omega));         % Preallocate k_0 array

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE));  % Preallocate
Cg = zeros(size(X_kRE));
a1 = zeros(size(X_kRE));

for ii = 1:length(omega)
    dispersion_eq = @(k) omega(ii)^2 - g * k * tanh(k * h);
     X_kRE(ii) = fzero(dispersion_eq, [0.01, 100]);  % ← WIDER INTERVAL FOR SCALED MODEL
end

for idx = 1:length(X_kRE)
    wk = X_kRE(idx);                % dimensional wavenumber  
    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);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %=============== 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);
    opts_lsq = optimoptions('lsqnonlin','Display','off');

    k = zeros(N,1);
    Znd = zeros(N,1);

    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;
            k(n) = lsqnonlin(fun_alpha, guess, x0_left, x0_right, opts_lsq);
            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 \lambda_m = m*pi/b
    m = zeros(M,1);
    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));

    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));
    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));
    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));
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%% Define Matrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    B=zeros(N,M);
    B(1,1) = (4*sinh(wk*b)) / (RE * wk*log(RE/RI) *cosh(wk*h));

    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

    for i=2:N
        B(i,1)=(4*sin(k(i)*b))/(RE*k(i)*log(RE/RI)*cos(k(i)*h));
    end

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

    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

    for i = 2:N
        C(i,1) = -4 * sin(k(i)*b) / (RE *k(i)* log(RE/RI) * cos(k(i)*h));
    end

    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

    % Matrix D
    D = zeros(M,N);
    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));

    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

    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

    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

    % Matrix E
    E = zeros(N, M);
    E(1,1) = (4 * sinh(wk*b)) / (RI *wk* log(RE/RI) * cosh(wk*h));

    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

    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

    % Matrix F
    F = zeros(N,M);
    F(1,1) = (-4 * sinh(wk*b)) / (RI * wk*log(RE/RI) * cosh(wk*h));

    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

    for i = 2:N
        F(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_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

    % Matrix W
    W = zeros(N,Q);
    for q = 1:Q
        r = exp(-2*chi(q)*b);
        coth_chib = (1 + r) / (1 - r);
        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 ) );
    end

    for n = 2:N
        for q = 1:Q
            r = exp(-2*chi(q)*b);
            coth_chib = (1 + r) / (1 - r);
            W(n,q) = -(4*chi(q)/cos(k(n)*h)) * ...
                      ( ( k(n)*sin(k(n)*b)*coth_chib + chi(q)*cos(k(n)*b) ) ...
                        / ( k(n)^2 + chi(q)^2 ) );
        end
    end

    % E1 and S_q
    f2 = omega(idx)^2 / g;
    E1  = zeros(Q,1);
    S_q = zeros(Q,1);

    for q = 1:Q
        exp2chi_d = exp(-2 * chi(q) * d);
        exp2chi_b = exp(-2 * chi(q) * b);
        exp2chi_h = exp(-2 * chi(q) * h);
        exp1chi_d = exp(-chi(q) * d);

        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);
        E1(q) = 2 * (BB * exp2chi_d - 1);

        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;
    end

    % Matrix H
    H1 = zeros(Q,1);
    for q = 1:Q
        r   = exp(-2*chi(q)*b);
        mid   = 4*r/(1 - r^2) - E1(q);
        H1(q) = mid + 1i*a1(idx)*chi(q)/omega(idx);
    end
    H = diag(H1);

    % Matrix G
    G = zeros(Q, N);
    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) );
    end

    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) );
        end
    end

    % Right hand side vector U
    U = zeros(2*M + 2*N + Q, 1);

    for i = 1:M
        if i == 1
            U(i, 1) = (sinh(wk*b))/((b*wk)*cosh(wk*h))*besselj(l, wk*RE) ;
        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);
        end
    end

    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);
        else
            U(j + M, 1) = 0;
        end
    end

    for i = 1:M
        U(i + M + N, 1) = 0;
    end

    for j = 1:N
        U(j + M + N + M, 1) = 0;
    end

    % Identity and zero matrices
    I_M = eye(M);
    I_N = eye(N);
    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 ];

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % NONLINEAR ITERATION LOOP
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    converged = false;
    psi = zeros(Q,1);
    T_old = [];

    for iter = 1:max_iter
        % (A) Update 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

        % (B) Solve the linear system
        T = S \ U;

        % (C) Extract coefficients
        b_vec = T(1:M);
        a_vec = T(M+1 : M+N);
        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);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % === MOLIN'S DISCRETIZATION FOR NONLINEAR INTEGRAL ===
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % STEP 1: Evaluate v_D at segment midpoints using CURRENT coefficients
        v_D_mid = zeros(1, Npt);
        for j = 1:Npt
            r = R_mid(j);
            % Call your v_D function at this radius
            v_D_mid(j) = v_D(N, Q, r, Z0d, wk, RI, l, d_vec, Znd, k, e_vec, chi);
        end

        % STEP 2: Compute nonlinear term F_j = |v_D| * v_D
        F_vec = abs(v_D_mid) .* v_D_mid;  % This handles complex values correctly

        % STEP 3: Compute psi(q) using analytical Bessel integration
        for q = 1:Q
            chi_q = chi(q);  % Current eigenvalue

            sum_integral = 0;
            for j = 1:Npt
                % Analytical integral: ∫[R_j to R_{j+1}] J_l(chi*r) * r dr
                % = (1/chi) * [R_{j+1}*J_{l+1}(chi*R_{j+1}) - R_j*J_{l+1}(chi*R_j)]

                if abs(chi_q) < 1e-10
                    % Special case: chi ≈ 0
                    integral_j = (R_seg(j+1)^2 - R_seg(j)^2) / 2;
                else
                    integral_j = (R_seg(j+1) * besselj(l+1, chi_q*R_seg(j+1)) - ...
                                 R_seg(j) * besselj(l+1, chi_q*R_seg(j))) / chi_q;
                end

                % Add contribution from this segment
                sum_integral = sum_integral + F_vec(j) * integral_j;
            end

            % Final psi value for this mode
            psi(q) = sum_integral;
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % === END OF MOLIN'S DISCRETIZATION ===
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % (E) Check convergence
        if ~isempty(T_old)
            diff_T = max(abs(T - T_old));
            if diff_T < tol
                converged = true;
                break
            end
        end
        T_old = T;

    end  % End iteration loop


    % Compute 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
fprintf("term1: %.3e, sum1: %.3e, phiUell: %.3e\n", term1, sum1, phiUell);
    eta0(idx) = abs(term1+sum1+phiUell);

end  % End frequency loop

term1: 4.535e-07, sum1: -1.221e-05, phiUell: -4.081e-07
term1: 1.599e-03, sum1: 2.875e-03, phiUell: 8.211e-04
term1: NaN, sum1: NaN, phiUell: NaN
term1: NaN, sum1: NaN, phiUell: NaN
term1: NaN, sum1: NaN, phiUell: NaN
term1: NaN, sum1: NaN, phiUell: NaN
term1: NaN, sum1: NaN, phiUell: NaN
term1: 9.423e-01, sum1: -1.384e-02, phiUell: -6.686e-03
term1: 9.526e-01, sum1: -1.403e-02, phiUell: -5.364e-03
term1: 9.603e-01, sum1: -1.307e-02, phiUell: -4.407e-03

% Plotting
fprintf('Length of T1: %d\n', length(T1));

Length of T1: 10

fprintf('Length of eta0: %d\n', length(eta0));

Length of eta0: 10

plot(T1, eta0, 'k', 'LineWidth', 1.5);

xlabel('$T$ (s)', 'Interpreter', 'latex');
ylabel('$|\eta / (iA)|$', 'Interpreter', 'latex');
title('Wave motion amplitude - 1:100 Scale Model', 'Interpreter', 'latex');
legend({'$\tau=0.2$ (Scale Model)'}, ...
       'Interpreter','latex','Location','northwest');

grid on;
xlim([0.5 3.5]);  % Scaled time range
ylim([0 4.5]);

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

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

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

    % Return the final value of v_D
    v_D_val = term1 + sum1 + sum2;
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 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 = 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 = 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 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

[Torsten]

The way you compute the integral is important for the runtime of your code, but is not responsible for and can't avoid the divergence. The fixed-point iteration for T in the part

 psi = zeros(Q,1);  
for iter = 1:max_iter
        % (A) Update 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
        % (B) Solve the linear system
        T = S \ U;
        % (C) Extract coefficients
        b_vec = T(1:M);
        a_vec = T(M+1 : M+N);
        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);
       % Somehow solve the integrals here
        psi = something
        % (E) Check convergence
        if ~isempty(T_old)
            diff_T = max(abs(T - T_old));
            if diff_T < tol
                converged = true;
                break
            end
        end
        T_old = T;
    end  % End iteration loop
    

diverges because the fixed-points for your second parameter case are non-attractive. This results in unrealistic values for d_vec = T(2*M+N+1:2*M+2*N) and e_vec = T(2*M+2*N+1:end) in the integration with the consequence that the integration also returns Inf or NaN values. Thus it won't help to modify the integration - it is ok with "integral" or "trapz" as is.

You will have to work on the method how to solve

S*T = U(T)

or

S*T-U(T) = 0

At the moment, you iterate T as

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

As I tried to explain in one of my comments from above (and as I proved with the code I use), the function

g(T) = S\U(T)

has eigenvalues of absolute value > 1 in the solution points T* of your system S*T - U(T) = 0 for certain parameter constellations. Thus when you use the fixed-point iteration for these parameters, you will never get T*.

Just to show you what I mean: For your code example from above, I get the following maximum absolute values of the eigenvalues of the iteration function in the solution point T* for the 10 values of X_kRE:

4.4e-3, 0.895, 6.8785, 19.356, 5.3709, 2.9084, 1.7406, 1.1185, 0.7607, 0.5418.

As you can see, from theory, you should only get convergence of the fixed-point method for the first and last two X_kRE values. You were lucky that the third last iteration converged, too.

Sometimes modifying the iteration function can solve such problems. E.g. if you want to solve

tan(x) = x

and iterate

x_(n+1) = tan(x_n)

you get divergence, but if you set x = atan(y) and write your equation as

y = atan(y)

and iterate

y_(n+1) = atan(y_n)

your fixed-point iteration converges (because the iteration function g(y) = atan(y) has derivative g'(y) = 1/(1+y^2) which is < 1).

Maybe you can find such a transformation for your case. Do you have an idea how exactly U(2*N+2*M+1:2*N+2*M+Q) look like as functions of T ? Are they quasi-quadratic in T ?

I didn't think about it more deeply and simply applied Newton's method to compute e_vec = T(2*M+2*N+1:end). For the other components of T you have linear equations and thus a_vec, b_vec, c_vec and d_vec can be deduced from e_vec.

I will have a look in the article later.

[Torsten]

I read the article and - correct me if I'm wrong - the problem considered in the article is instationary while yours is stationary. Thus the numerical tool used in the article is a solver for ordinary differential equations while your tool must be a solver for algebraic equations. I don't see how to extract useful information for your case.

[Javeria]

@Torsten the only difference is that they have time dependent study and i have the time independent they also used the same method and form the system of equations and for the non linear integral they mentioned that method.

[Torsten]

the only difference is that they have time dependent study and i have the time independent they also used the same method and form the system of equations and for the non linear integral they mentioned that method.

I know this. The authors of the article will most probably use a scientific code that solves their system of ordinary differential equations. In MATLAB, this would be "ode15s". For your problem, MATLAB offers "fsolve" to solve nonlinear systems of algebraic equations. It's unusual that authors develop their own solution method (like you try to do) - they use general-purpose software because it's mostly faster and more reliable. And concerning the integral: you can use the author's method if you like, but it won't change anything in the behaviour of your fixed-point iteration.

Let me put it clear:

Newton's method is the standard method to solve systems of nonlinear equations. If you find an article that solves some kind of stationary problem like you do, you will find with 99% probability that the authors used Newton's method. And it works as you can see from the code I use. If you only want to use it for your own work and don't want to sell it, speed doesn't matter in my opinion as long as you get results in an acceptable time. For small values for Q, you get results almost instantly. If you run it for larger values, you can start several jobs in parallel and get the results next day. For me, this is acceptable. If you further want to speed up the code, start anew and use Fortran or C.

If you have questions concerning the code or an idea how to use another method instead of your fixed-point iteration, you are again welcome to ask.

[Javeria]

I find an article where they solved the non linear pressure drop condition. Yes offcouurse the expression are different but overall if we speak generally then the method to solve this kind of problem should be same, below i have another idea from an article to solve this kind of problem. They mentioned the 6 points guassian quadrature integration. and then using the newton raphson method.

[Torsten]

The code I provided uses Newton-Raphson to solve for e_vec and Gauss-Kronos quadrature to determine the integral. Why do you think Newton-Raphson for (24) (whatever it is) and 6-point Gaussian quadrature will be better ?

If you find references that use Newton-Raphson, the respective code will behave similarly in speed as the one I provided.

[Paul]

@Torsten

Referring to this comment where it's stated

"the function

g(T) = S\U(T)

has eigenvalues of absolute value > 1"

If the iteration function is

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

then shouldn't we be concerned about eigenvlaues of

inv(S)

rather than

g(T) = S\U(T) ?

I'm sure that's how it should be if U(T) were simply U(T) = T.

[Torsten]

Sorry, I meant the eigenvalues of the Jacobian of the iteration function

g(T) = S\U(T)

with respect to T here.

If U(T) = T, this will be

inv(S)

If U is a general function of T, it should be

inv(S) * Jac(U)

I guess.

[Javeria]

@Torsten can you please share that code which only use the newton raphson and the Gauss-Kronos quadrature to determine the integral. Thanks

[Torsten]

Newton-Raphson is used for solution_method = 1, 2, 3 or 4.

Gauss-Kronos quadrature is used for integral_method = 1.

Solver_method = 12, 22 or 32 can be selected.

I recommend solution_method = 2, integral_method = 2 (the trapoidal rule gave almost identical results to Gauss-Kronos, but was a lot faster for bigger problems) and solver_method = 32 (nleq1 was the fastest of the Newton-solvers).

I don't want to destroy choices because I'm a pessimist and want to have a plan B if plan A doesn't succeed.

The abscissa for plotting is computed as RE * X_kRE * tanh(X_kRE * h). I saw that you changed this to

2*pi / sqrt(g * X_kRE * tanh(X_kRE * h)) in your last posted code.

[Torsten]

@Javeria

In the meantime, I made slight modifications to the code. If you plan to use it, start from the attached version.

[Javeria]

@Torsten I successfully run your code thanks but the plot i have this blue line also. What does it show?

[Torsten]

The blue curve is the curve that is created at the X_kRE values you prescribe.

The solution is advanced from one X_kRE value to the next, and the initial guess for the vector of unknowns at step i is taken as the solution at step (i-1). Sometimes it may happen that the solver doesn't converge at step i. Then a value for your X_kRE is inserted in the middle of the interval [X_kRE(i-1),X_kRE(i)] and the solver is called for (X_kRE(i-1)+X_kRE(i))/2, again with the solution at step (i-1) as initial guess. The red curve shows the solution in your prescribed values for X_kRE plus the solution in the points for X_kRE that were inserted.

I write two files to your working directory: data_test and data_test_full. data_test contains the solution in your prescribed X_kRE vector, data_test_full contains the solution in your prescribed X_kRE vector plus in the points that were inserted. You can use these files later to make plots for different cases. Don't forget to rename the files after a run has finished - otherwise they will be overwritten in the next run.

Further there is an option to restart the computation from the results obtained at a certain value for X_kRE if a run takes too long and you want to shut down your computer. The name of the restart file is solution_test.

You can prescribe the names of all the files written in the section "Solution file names" of the code.

[Javeria]

@Torsten Thanks for your valuable suggestions. Below is the script to solve the problem and i did validation . I tested the below code for different N upto 100 it converges and the results are similar with the full scale model also.

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 = 34.5;                 % Inner radius (m)
% % ratio = 47.5/34.5;         % Ratio R_E / R_I
% % % ratio = 2.0;
% % RE = ratio*RI;             % outer radius (m)
% % h = 200/34.5 * RI;         % Water depth (m)
% % d = 38/34.5 * RI;              % Interface depth (m)
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% === SCALED GEOMETRIC PARAMETERS (divided by 100) ===
RI = 0.345;                % Inner radius (m) [34.5/100]
ratio = 47.5/34.5;         % Ratio R_E / R_I (dimensionless - unchanged)
RE = ratio*RI;             % outer radius (m) [0.4750 m]
h = 200/34.5 * RI;         % Water depth (m) [2.0 m]
d = 38/34.5 * RI;          % Interface depth (m) [0.38 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
%b1 =0;
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% T1 = linspace(5, 35, 10);           % Define wave period range [s]
T1 = linspace(0.5, 3.5, 50);        % Scaled period range [s] [5/10 to 35/10]
fprintf('Check T1 immediately after linspace: %d\n', length(T1));

Check T1 immediately after linspace: 50

omega = 2 * pi ./ T1;               % Angular frequency [rad/s]
X_kRE = zeros(size(omega));        % Preallocate k_0 array
 % % X_kRE = linspace(0.0025, 0.16,100);         %%%%%% this is now k_0 
% kh_range = linspace(0.01, 10, 100);   % Covering shallow to deep regimes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ==============Initialize eta0 array before loop
%===============Initialize omega array before loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta0 = zeros(size(X_kRE));  % Preallocate
for ii = 1:length(omega)
    dispersion_eq = @(k) omega(ii)^2 - g * k * tanh(k * h);
    X_kRE(ii) = fzero(dispersion_eq, [0.001, 100]);
end
for idx = 1:length(X_kRE)
    % if I_given_wavenumber == 1
        wk = X_kRE(idx);                % dimensional wavenumber  
        %omega(idx) = sqrt(g * wk * tanh(wk * h)); % dispersion relation
       % fun_alpha = @(x) omega(idx)^2 + g*x.*tan(x*h);     %dispersion equation for evanescent modes

% % 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 (omega=%6.3f s) ---\n', idx, 2*pi/T1(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 ];
fprintf('cond(S) = %.2e\n', cond(S));

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               NEWTON–RAPHSON NONLINEAR SOLVER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% --- 6-point Gauss–Legendre quadrature on [0, RI]
xg = [-0.9324695142; -0.6612093865; -0.2386191861; ...
       0.2386191861;  0.6612093865;  0.9324695142];
wg = [ 0.1713244924;  0.3607615730;  0.4679139346; ...
       0.4679139346;  0.3607615730;  0.1713244924];

r_g = 0.5 * RI * (xg + 1);
w_g = 0.5 * RI * wg;

% --- Newton parameters
psi = zeros(Q,1);            % initial guess (linear solution)
tol_NR = 1e-6;
max_NR = 100;

for it = 1:max_NR

    % ===============================
    % 1) Evaluate F(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

    % Solve linear system
    T = S \ U;

    d_vec = T(2*M+N+1 : 2*M+2*N);
    e_vec = T(2*M+2*N+1 : end);

    % Compute nonlinear residual F
    Fnl = zeros(Q,1);
    for q = 1:Q
        Iq = 0;
        for kq = 1:6
            r = r_g(kq);
            v = v_D(N,Q,r,Z0d,wk,RI,l,d_vec,Znd,k,e_vec,chi);
            Iq = Iq + w_g(kq) * abs(v) * v ...
                       * besselj(l,chi(q)*r) * r;
        end
        Fnl(q) = psi(q) - Iq;
    end

    % Convergence check
    if norm(Fnl,inf) < tol_NR
        fprintf('Newton converged in %d iterations\n', it);
        break
    end

    % ===============================
    % 2) Jacobian (finite difference)
    % ===============================
    J = zeros(Q,Q);
    eps_fd = 1e-6;

    for j = 1:Q
        psi_p = psi;
        psi_p(j) = psi_p(j) + eps_fd;

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

        T_p = S \ U;

        d_p = T_p(2*M+N+1 : 2*M+2*N);
        e_p = T_p(2*M+2*N+1 : end);

        for q = 1:Q
            Iqp = 0;
            for kq = 1:6
                r = r_g(kq);
                v = v_D(N,Q,r,Z0d,wk,RI,l,d_p,Znd,k,e_p,chi);
                Iqp = Iqp + w_g(kq) * abs(v) * v ...
                            * besselj(l,chi(q)*r) * r;
            end
            J(q,j) = (psi_p(q) - Iqp - Fnl(q)) / eps_fd;
        end
    end

    % ===============================
    % 3) Newton update
    % ===============================
    delta = J \ Fnl;
    psi = psi - delta;

end

if it == max_NR
    warning('Newton did not converge');
end





% % ---- 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

--- idx   1 (omega=12.566 s) ---

cond(S) = 7.06e+04

Newton converged in 1 iterations

--- idx   2 (omega=11.195 s) ---

cond(S) = 4.73e+04

Newton converged in 1 iterations

--- idx   3 (omega=10.094 s) ---

cond(S) = 3.45e+04

Newton converged in 1 iterations

--- idx   4 (omega= 9.190 s) ---

cond(S) = 2.64e+04

Newton converged in 1 iterations

--- idx   5 (omega= 8.435 s) ---

cond(S) = 2.10e+04

Newton converged in 2 iterations

--- idx   6 (omega= 7.794 s) ---

cond(S) = 1.71e+04

Newton converged in 3 iterations

--- idx   7 (omega= 7.244 s) ---

cond(S) = 1.42e+04

Newton converged in 4 iterations

--- idx   8 (omega= 6.767 s) ---

cond(S) = 1.20e+04

Newton converged in 4 iterations

--- idx   9 (omega= 6.348 s) ---

cond(S) = 1.02e+04

Newton converged in 5 iterations

--- idx  10 (omega= 5.978 s) ---

cond(S) = 8.82e+03

Newton converged in 9 iterations

--- idx  11 (omega= 5.649 s) ---

cond(S) = 7.68e+03

Newton converged in 14 iterations

--- idx  12 (omega= 5.354 s) ---

cond(S) = 6.76e+03

Newton converged in 19 iterations

--- idx  13 (omega= 5.089 s) ---

cond(S) = 6.02e+03

Newton converged in 27 iterations

--- idx  14 (omega= 4.848 s) ---

cond(S) = 5.56e+03

Newton converged in 30 iterations

--- idx  15 (omega= 4.630 s) ---

cond(S) = 5.95e+03

Newton converged in 30 iterations

--- idx  16 (omega= 4.430 s) ---

cond(S) = 8.70e+03

Newton converged in 11 iterations

--- idx  17 (omega= 4.247 s) ---

cond(S) = 1.67e+04

Newton converged in 10 iterations

--- idx  18 (omega= 4.078 s) ---

cond(S) = 5.87e+03

Newton converged in 26 iterations

--- idx  19 (omega= 3.922 s) ---

cond(S) = 3.26e+03

Newton converged in 39 iterations

--- idx  20 (omega= 3.778 s) ---

cond(S) = 2.97e+03

Newton converged in 36 iterations

--- idx  21 (omega= 3.644 s) ---

cond(S) = 2.78e+03

Newton converged in 24 iterations

--- idx  22 (omega= 3.519 s) ---

cond(S) = 2.62e+03

Newton converged in 20 iterations

--- idx  23 (omega= 3.402 s) ---

cond(S) = 2.49e+03

Newton converged in 19 iterations

--- idx  24 (omega= 3.293 s) ---

cond(S) = 2.39e+03

Newton converged in 19 iterations

--- idx  25 (omega= 3.190 s) ---

cond(S) = 2.31e+03

Newton converged in 18 iterations

--- idx  26 (omega= 3.094 s) ---

cond(S) = 2.25e+03

Newton converged in 17 iterations

--- idx  27 (omega= 3.004 s) ---

cond(S) = 2.21e+03

Newton converged in 16 iterations

--- idx  28 (omega= 2.918 s) ---

cond(S) = 2.18e+03

Newton converged in 15 iterations

--- idx  29 (omega= 2.838 s) ---

cond(S) = 2.16e+03

Newton converged in 14 iterations

--- idx  30 (omega= 2.761 s) ---

cond(S) = 2.15e+03

Newton converged in 13 iterations

--- idx  31 (omega= 2.689 s) ---

cond(S) = 2.15e+03

Newton converged in 12 iterations

--- idx  32 (omega= 2.620 s) ---

cond(S) = 2.16e+03

Newton converged in 11 iterations

--- idx  33 (omega= 2.555 s) ---

cond(S) = 2.17e+03

Newton converged in 11 iterations

--- idx  34 (omega= 2.493 s) ---

cond(S) = 2.18e+03

Newton converged in 10 iterations

--- idx  35 (omega= 2.434 s) ---

cond(S) = 2.20e+03

Newton converged in 9 iterations

--- idx  36 (omega= 2.377 s) ---

cond(S) = 2.22e+03

Newton converged in 9 iterations

--- idx  37 (omega= 2.324 s) ---

cond(S) = 2.25e+03

Newton converged in 8 iterations

--- idx  38 (omega= 2.272 s) ---

cond(S) = 2.27e+03

Newton converged in 8 iterations

--- idx  39 (omega= 2.223 s) ---

cond(S) = 2.30e+03

Newton converged in 7 iterations

--- idx  40 (omega= 2.176 s) ---

cond(S) = 2.37e+03

Newton converged in 7 iterations

--- idx  41 (omega= 2.131 s) ---

cond(S) = 2.46e+03

Newton converged in 7 iterations

--- idx  42 (omega= 2.087 s) ---

cond(S) = 2.55e+03

Newton converged in 6 iterations

--- idx  43 (omega= 2.046 s) ---

cond(S) = 2.65e+03

Newton converged in 6 iterations

--- idx  44 (omega= 2.006 s) ---

cond(S) = 2.74e+03

Newton converged in 6 iterations

--- idx  45 (omega= 1.967 s) ---

cond(S) = 2.84e+03

Newton converged in 6 iterations

--- idx  46 (omega= 1.930 s) ---

cond(S) = 2.94e+03

Newton converged in 5 iterations

--- idx  47 (omega= 1.895 s) ---

cond(S) = 3.03e+03

Newton converged in 5 iterations

--- idx  48 (omega= 1.860 s) ---

cond(S) = 3.13e+03

Newton converged in 5 iterations

--- idx  49 (omega= 1.827 s) ---

cond(S) = 3.23e+03

Newton converged in 5 iterations

--- idx  50 (omega= 1.795 s) ---

cond(S) = 3.32e+03

Newton converged in 5 iterations

% 
% % % ---------- Summary of iterations for T in [5, 35] s ----------
% % T_all = 2*pi./omega(:);
% % mask  = (T_all >= 5) & (T_all <= 35);
% % [Ts, order] = sort(T_all(mask), 'descend');   % sort by T (optional)
% % idx_list    = find(mask);
% % idx_list    = idx_list(order);
% % iters_list  = iters_used(mask);
% % iters_list  = iters_list(order);
% 
% % fprintf('\nIterations for T in [5, 35] s (sorted by T):\n');
% % fprintf('  idx   T [s]   iters\n');
% % for i = 1:numel(idx_list)
% %     fprintf(' %4d  %7.3f     %2d\n', idx_list(i), Ts(i), iters_list(i));
% % end
% 
% % % (optional) write the same table to a CSV file
% % tbl = table(idx_list, Ts, iters_list, eta0(idx_list)', ...
% %     'VariableNames', {'idx','T','iters','eta0'});
% % writetable(tbl, 'iters_5to35.csv');
% % fprintf('Wrote iters_5to35.csv with %d rows.\n', height(tbl));
% 
% % % Show which index is the resonance and how many iterations it used
% % [~, i_peak] = max(eta0);
% % fprintf('\nRES peak at idx %d (T=%6.3f s): iters=%d\n', ...
% %         i_peak, 2*pi/omega(i_peak), iters_used(i_peak));
% % 
% % % (optional) print a small window around resonance
% % win = 5;  % ±5 indices around the peak
% % i0  = max(1, i_peak - win);
% % i1  = min(length(X_kRE), i_peak + win);
% % fprintf('\nIterations near resonance (idx %d):\n', i_peak);
% % fprintf('  idx   T [s]   iters\n');
% % for j = i0:i1
% %     fprintf(' %4d  %7.3f     %2d\n', j, 2*pi/omega(j), iters_used(j));
% % end
% 
% % % (optional) quick plot: iterations vs T
% % figure(1); clf;
% % plot(T_all, iters_used, '-o'); grid on;
% % xlabel('T [s]'); ylabel('iterations'); title('Nonlinear iterations vs T');
% % xlim([5 35]);
% 
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % %                   Load CSV files
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % Load CSV files HERE (after the loop ends)
%     % % data_chen = load('dr_chen_no_dissp.csv');
% % % data_liu  = load('dr_liu.csv');
%      % data_exp  = load('experimental.csv');
% % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % % %              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 
fprintf('Length of T1: %d\n', length(T1));

Length of T1: 50

fprintf('Length of eta0: %d\n', length(eta0));

Length of eta0: 50

plot(T1, eta0, 'k', 'LineWidth', 1.5); %%% 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');

Warning: Ignoring extra legend entries.

% 
  grid on;
  % xlim([5 35]);  % Match the reference plot
  xlim([0.5 3.5]);  % Scaled time range
  ylim([0 4.5]);  % Optional: based on expected peak height

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

Time consuming =  21.406 s

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

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

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

    % Return the final value of v_D
    v_D_val = term1 + sum1 + sum2;
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 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 = 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 = 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 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