how to evaluate the integral of the expression involving bessel functions.

I have the following expression which i need to find.

then how i can find the integration? here this k is vector which contain k0,k1,k2, and same d and e is also a vector which contain d1,d2,d3 and e1, e2 ,e3 . and is also given. Then how i can evaluate this integration.

[Torsten]

Why is psi a function of r ?

Which index p is meant when you write nu_D(r)*abs(nu_D(r))*besselj(chi_p^l * r) r as the integrand ?

Do you have functions to compute the derivatives of besselj and besseli ?

Why does your vector k has only 3 elements instead of 4 for k0,k1,k2 and k3 ?

[Javeria]

@Torsten yes i know the functions that how we can compute derivatives of bessel functions here we have k0,k1,k2 belongs to N=3.

Answer 1

[Shishir Reddy]

Hi @Javeria

Kindly refer the following steps to numerically evaluate the given integral using MATLAB.

1. Define the Parameters:

% Given values
RI = 0.21;
k = [0.281, 0.2, 0.1];
d = [0.1, 0.3, 0.4];   % d0, d1, d2, d3 → d(1), d(2), d(3), d(4)
e = [0.3, 0.09, 0.7];
chi = [0.01, 0.03, 0.02];

2. Define ν_D(r):

nu_D = @(r) ...
    d(1) * besselj(0, k(1)*r) / besselj(0, k(1)*RI) + ...
    d(2) * besseli(1, k(2)*r) / besseli(1, k(2)*RI) + ...
    d(3) * besseli(1, k(3)*r) / besseli(1, k(3)*RI) + ...
    e(1) * besselj(1, chi(1)*r) / besselj(1, chi(1)*RI) + ...
    e(2) * besselj(1, chi(2)*r) / besselj(1, chi(2)*RI) + ...
    e(3) * besselj(1, chi(3)*r) / besselj(1, chi(3)*RI);

3. Define the Integrand Function:

chi1 = chi(1);
integrand = @(r) abs(nu_D(r)).^2 .* besselj(0, chi1 * r);

4. Use 'integral' to compute the result:

psi = integral(integrand, 0, RI);
disp(['ψ = ', num2str(psi)]);

For more information regarding 'integral' function, kindly refer the following documentation -

https://www.mathworks.com/help/matlab/ref/integral.html

I hope this helps.

[Javeria]

@Shishir Reddy hello thanks for your response i might ask this question that how we can write this v_D(r) the way you write i understand but if we have a large N. then it's very difficult to write it.

Answer 2

[Torsten]

l = 0;
R_I = 2;
N = 3;
P = 3;
k = [0.28,0.2,0.1,0.05];   %[k0,k1,k2,k3]
e = [0.3,0.9,0.7];         %[e1,e2,e3]
d = [0.1,0.3,0.4,0.5];     %[d0,d1,d2,d3]
chi = [0.01,0.03,0.02];    %[chi1,chi2,chi3]
format long
psi = integral(@(r)fun(r,l,R_I,N,P,k,e,d,chi),0,R_I)

psi = 
     3.506800645644731e+02

function values = fun(r,l,R_I,N,P,k,e,d,chi)
  nu_D_values = nu_D(r,l,R_I,N,P,k,e,d,chi);
  chip = chi(end);          %-> which p-index from p has to be taken here ?
  values = nu_D_values.*abs(nu_D_values).*besselj(l,chip*r).*r;   
end

function nu_D_values = nu_D(r,l,R_I,N,P,k,e,d,chi)
  nu_D_values = d(1)*besselj(l,k(1)*r)/dbesselj(l,k(1)*R_I);
  for n = 1:N
    nu_D_values = nu_D_values + d(n+1)*besseli(l,k(n+1)*r)/dbesseli(l,k(n+1)*R_I);  
  end
  for p = 1:P
    nu_D_values = nu_D_values + e(p)*chi(p)*besselj(l,chi(p)*r)/dbesselj(l,chi(p)*R_I);
  end
end

function value = dbesseli(nu,z)
   value = 0.5*(besseli(nu-1,z) + besseli(nu+1,z));
   %value = nu./z.*besseli(nu,z) + besseli(nu+1,z);
end

function value = dbesselj(nu,z)
   value = 0.5*(besselj(nu-1,z) - besselj(nu+1,z));
   %value = nu./z.*besselj(nu,z) - besselj(nu+1,z);
end

[Javeria]

@Torsten Thanks. However i having a confusion regarding this v_D(r) as we have summation for the 2nd and third term so can we define it like

%%%% for the 2nd term
sum1 =0;
for n=2:N
sum1 = sum1+d(n)*besseli(l,k(n)*r)/dbesseli(l,k(n)*R_I);
end

For the summation term should i defined my terms in this way or the above one.

[Javeria]

@Shishir Reddy@Torsten I write in this way

function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, dn, Znd, k, eq, chi)
% Computes v_D(r) for arbitrary r, given coefficients and parameters.
% Inputs:
%   N  - number of n modes
%   Q  - number of q modes
%   r  - radial coordinate (scalar)
%   Z0d, wk, RI, l, dn, Znd, k, eq, chi - 

    % First term: n=0 mode
   term1 = (dn(1) * Z0d * besselj(l, wk*r)) / (wk*dbesselj(l, wk*RI));

    % Second term: sum over n=2:N
    sum1=0;
    for nidx = 2:N % For N=3, this runs for nidx = 2 and nidx = 3
        sum1= sum1 + dn(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) /(k(nidx)* dbesseli(l, k(nidx)*RI));
    end

    % Third term: sum over q=1:Q
    sum2 =0;
    for qidx = 1:Q
        sum2= sum2+ eq(qidx) * chi(qidx) * besselj(l, chi(qidx)*r) / (chi(qidx)*dbesselj(l, chi(qidx)*RI));
    end
 % Sum all terms to get the final v_D(r) value
    v_D_val = term1 + sum1 + sum2;
end

and then calling it

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

[Torsten]

If you have a function f(x) and you want to take the derivative of f at k*x (thus f'(k*x)), you just evaluate f'(x) at k*x.

If you have a function g(x) = f(k*x) and you want to take the derivative of g at x, it's necessary to use the chain rule - thus g'(x) = k*f'(k*x).

In your case, you want to take the derivative of f at k*x, not the derivative of g at x. Thus the inner derivatives (wk, k and chi) in the denominator of your expressions are wrong in my opinion.

You can also make the dimensional test to see which version is correct: wk*dbesselj(l, wk*RI), k(nidx)* dbesseli(l, k(nidx)*RI) and chi(qidx)*dbesselj(l, chi(qidx)*RI) have unit 1/m, dbesselj(l, wk*RI), dbesseli(l, k(nidx)*RI) and dbesselj(l, chi(qidx)*RI) have unit 1.

And I would handle [Z0d,Znd] and [wk,k] the same way as you handle d. In your code, you don't reference the first elements of Znd and k, but you could put Z0d and wk at this first position and remove them from the input list. Or you could define d0 and d and additionally put d0 in the input list. The main thing is: handle the cases consistently.

And I wouldn't put it as

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

because this way, your function v_D has to be evaluated two times. Use "integrand" as a function as I did, not as a function handle. The index q for which you want to evaluate chi can easily be passed to this function.

But all in all, you should arrange the code in a way that you understand it - you are the one who has to debug it at the final stage. Especially there is one thing you should remember: use the same variable names in your documentation and in your code - it will make it much easier to compare what you coded and what you should have coded.

[Javeria]

@Torsten Thanks for your response, I can share my full code. Actually i want to solve my system of equations for unknown coefficients b0,bm,a0,an,c0,cm,d0,dn,eq. However total i have 9 equations out of which 1 have non linear term due to this psi term. So i treat my this non linear term iteratively to solve my system of equations. However when i run my code i get zero vector for my unknowns which is physically not correct.

In my code i want to treat my psi iteratively that in 1st iteration it take the initial guess which i set to zero. And then onward it is updated by this expression. Below is my main script where i am calling helper functions for derivative of bessel functions and roots of bessel functions i set l=0 for the bessel functions mode.

The way i did for this psi is correct or not in my code?

% --- Main Script: Solve Non-Linear Moonpool Problem Iteratively ---
% Clear workspace, close all figures, and clear command window
close all;
clear all;
clc;
% --- 1. Define Global Constants and Parameters ---
% These define the truncation limits for the series summations in the equations.
N = 4;                   % Number of terms for 'n' modes (e.g., in Region 1/sum for 'a_n')
M = 4;                   % Number of terms for 'm' modes (e.g., in Region 2/sum for 'b_m')
Q = 3;                   % Number of terms for 'q' modes (e.g., in Region 3/sum for 'e_q')
% --- Physical Parameters ---
% These values are taken from the reference "Sphaier et al. (2007)".
RE = 47.5; % External radius of the main hull (m)
RI = 34.5; % Internal radius (m) - openness at the bottom of the moonpool
h = 200; % Water depth (m)
d = 38; % Draft (m)
g = 9.81; % Acceleration due to gravity (m/s^2)
b = (h-d); %  'b' as (h-d)
% --- Additional Constants ---
% These are constants specific to problem formulation.
gamma = 1.0;             % empirical discharge coefficient term is in the range [0.5, 1.0]
tau = 0.2;               % porosity ratio

% --- Bessel Function Order ---
% 'l' is the order of the Bessel functions used throughout your system of equations.
l = 0; % You can change this value as desired/required by the problem.

% This constant is derived from the formula (1 - tau) / (2 * gamma * tau^2)
b1 = (1 - tau) / (2 * gamma * tau^2);
% --- Iterative Solver Parameters ---
max_iter = 100;
tol = 1e-4;
converged = false;

% --- Wave Number and Group Velocity Calculation ---
X_kRE = linspace(5, 35, 100);   % X-axis (kRE) matching the plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I_given_wavenumber = 1;         % Example flag, as in your code
% Loop through each dimensionless wavenumber (k0*RE) value
% For each X_kRE(i), we calculate the corresponding k0 (wk), omega, Cg, and
% the evanescent modes k(n)
for i = 1:length(X_kRE)
% Calculate the progressive wave mode wavenumber (k0),
     %if I_given_wavenumber == 1
      % wk = k0 = (k0*RE) / RE
        wk    = X_kRE(i) / RE;                  % wavenumber
       omega = sqrt(g * wk * tanh(wk * h));     % wave frequency
Cg = (g*tanh(wk*h) + g*wk*h*(sech(wk*h))^2)*omega/(2*g*wk*tanh(wk*h));  %group velocity
        a1 = 0.93 * (1 - tau) * Cg;  % 
fun_alpha = @(x) omega^2 + g * x * tan(x * h); % dispersion eqn
   % end

    for n = 1:N
        if n == 1
            k(n) = -1i * wk;
        else
            x0_n = [(2*n - 3)*pi/2/h, (2*n - 1)*pi/2/h];
            k(n) = fzero(fun_alpha, mean(x0_n));
        end
    end
%%%%%%%%%%%%% Find derivative of Z0 at z=-d%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Z0d = (cosh(wk*h) * wk * sinh(wk*b))/(2*wk*h + sinh(2*wk*h));
% 
%%%%%%%%%%%%% Find derivative of Zn at z=-d%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Znd = zeros(1, N); % or however many n you have
for n = 2:N   % Typically n starts from 2 for these modes
    Znd(n) = -k(n) * (cos(k(n)*h) * sin(k(n)*b)) / (2*k(n)*h + sin(2*k(n)*h));
end

for m=2:M
    lambda(m) = ((m-1) * pi) / b;
end
% % Predefine identity and zero submatrices
O1  = zeros(M, M);  % O1 == O6, size M×M
O2  = zeros(M, N);  % O2 == O7, size M×N
O3  = zeros(M, Q);  % O3 == O8 size M×Q
O4  = zeros(N, N);  % O4 == O9, size N×N
O5  = zeros(N, Q);  % size N×Q
O10  = zeros(Q, M);  % O10 == O12, size Q×M
O11 = O5';          % O11 is the transpose of O5
I_M = eye(M);    % size M×M
I_N = eye(N);    % size N×N

%%%%%%%%%%%%%%%%%%%% Set Matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    A = zeros(M,N);
    % For first element A00
A(1,1) = (cosh(wk*h) * sinh(wk*b)) / (wk*b * (2.*wk*h + sinh(2.*wk*h))) ...
         * (besselh(l, 1, wk*RE) / wk*dbesselh(l, wk*RE));
% For 2nd Element A0n
for n = 2:N
    A(1,n) =  (cos(k(n)*h) * sin(k(n)*b)) / (k(n)*b * (2.*k(n)*h + sin(2.*k(n)*h)))...
        * (besselk(l, k(n)*RE) /k(n)* dbesselk(l, k(n)*RE));

end
% for element Am0
for m = 2:M
  
A(m,1) = 2. * (cosh(wk*h) * (-1)^(m-1) * wk * sinh(wk*b)) / b*((wk^2 + (lambda(m))^2) ...
    * (2.*wk*h + sinh(2.*wk*h))) * (besselh(l, 1, wk*RE) / wk*dbesselh(l, wk*RE));
end
%%%%%%%% for element Amn
for m = 2:M
    for n = 2:N
   
A(m,n) = 2.* (cos(k(n)*h) * (-1)^(m-1)* k(n) * sin(k(n)*b)) / b*((k(n)^2 - (lambda(m))^2)...
    * (2.*k(n)*h + sin(2.*k(n)*h))) * (besselk(l, k(n)*RE) / k(n)*dbesselk(l, k(n)*RE));

    end
end
%%%%%%%%%%%%%%%%%% Matrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 B = zeros(N, M); % Or whatever size you require
%%%%% Define B0%%%%%
B(1,1) = (4. * sinh(wk * b)) / (RE * log(RE / RI) * cosh(wk* h));
%%%%%%%% define Bom%%%%
for m = 2:M
    Rml_prime_RE = lambda(m) * (besselk(l, lambda(m)*RI) * dbesseli(l, lambda(m)*RE) ...
  - besseli(l, lambda(m)*RI) * dbesselk(l, lambda(m)*RE))...
  /(besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
    B(1,m) = Rml_prime_RE * (4 * wk^2 * (-1)^(m-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + (lambda(m))^2));
end
%%%% define Bn0%%%%
for n = 2:N
 B(n,1) = (4* sin(k(n)*b)) / (RE * log(RE/RI) * cos(k(n)*h));
end
%%%%%% Define Bnm %%%%%%%%
for n = 2:N
    for m = 2:M
        Rml_prime_RE = lambda(m) * (besselk(l, lambda(m)*RI) * dbesseli(l, lambda(m)*RE) ...
  - besseli(l, lambda(m)*RI) * dbesselk(l, lambda(m)*RE))...
  /(besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
        B(n,m) = Rml_prime_RE * (4* (k(n))^2 * (-1)^(m-1) * sin(k(n)*b)) / (cos(k(n)*h) * ((k(n))^2 - (lambda(m))^2));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix C %%%%%%%%%%%%%%%%%%%%%%%%
C = zeros(N,M);
C(1,1) = -4 * sinh(wk*b) / (RE * log(RE/RI) * cosh(wk*h));
% %%%%% DEfine C0m%%%%%%
for m = 2:M
     Rml_prime_star_RE = lambda(m) * (besseli(l, lambda(m)*RE) * dbesselk(l, lambda(m)*RE) ...
  - besselk(l, lambda(m)*RE) * dbesseli(l, lambda(m)*RE))...
  /(besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
    C(1,m) = Rml_prime_star_RE * (4 * wk^2 * (-1)^(m-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + (lambda(m))^2));
end
%%%%%%% Define Cn0%%%%%%
for n = 2:N
    C(n,1) = -4 * sin(k(n)*b) / (RE * log(RE/RI) * cos(k(n)*h));
end
%%%%%% Define Cnm%%%%%%
for n = 2:N
    for m =2:M

     Rml_prime_star_RE = lambda(m) * (besseli(l, lambda(m)*RE) * dbesselk(l, lambda(m)*RE) ...
  - besselk(l, lambda(m)*RE) * dbesseli(l, lambda(m)*RE))...
  /(besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
    C(n,m) = Rml_prime_star_RE * (4 * (k(n))^2 * (-1)^(m-1) * sin(k(n)*b)) / (cos(k(n)*h) * ((k(n))^2 - (lambda(m))^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) / wk*dbesselj(l, wk*RI));
%%%%% write D0n %%%%%%%%%%
for n =2:N
    D(1,n) = (cos(k(n)*h) * sin(k(n)*b)) / (k(n)*b * (2.*k(n)*h + sin(2.*k(n)*h))) * (besseli(l, k(n)*RI) / k(n)*dbesseli(l, k(n)*RI));

end
%%%%%% write Dm0 %%%%%%%%%
for m = 2:M
    D(m,1) = (2 * cosh(wk*h) * (-1)^(m-1) * wk * sinh(wk*b)) /(b * (2.*wk*h + sinh(2.*wk*h)) * (wk^2 + (lambda(m))^2))...
        *(besselj(l, wk*RI) /wk* dbesselj(l, wk*RI));
end
%%%%% Define Dmn%%%%%%%%
for m = 2:M
    for n = 2:N
        D(m,n) = (2 * cos(k(n)*h) * (-1)^(m-1) * k(n) * sin(k(n)*b)) /(b * (2.*k(n)*h + sin(2.*k(n)*h)) * ((k(n))^2 - (lambda(m))^2))...
        *(besseli(l, k(n)*RI) /k(n)* dbesseli(l, k(n)*RI));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%% Define Matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = zeros(N, M);  % Preallocate the matrix E
%%%%% Define E00%%%%%%%%%%%%
E(1,1) = (4 * sinh(wk*b)) / (RI * log(RE/RI) * cosh(wk*h));
%%%% Define Eom %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for m = 2:M
    Rml_prime_RI = lambda(m)*(besselk(l, lambda(m)*RI) * dbesseli(l, lambda(m)*RI) ...
  - besseli(l, lambda(m)*RI) * dbesselk(l, lambda(m)*RI))...
  / (besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
E(1,m) = Rml_prime_RI * (4 * wk^2 * (-1)^(m-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + lambda(m)^2));

end
%%%%%% Define Eno%%%%%%%%%%%%%%
for n = 2:N
    E(n,1) = (4 * sin(k(n)*b)) / (RI * log(RE/RI) * cos(k(n)*h));
end
for n = 2:N
    for m = 2:M
        Rml_prime_RI = lambda(m)*(besselk(l, lambda(m)*RI) * dbesseli(l, lambda(m)*RI) ...
  - besseli(l, lambda(m)*RI) * dbesselk(l, lambda(m)*RI))...
  / (besselk(l, lambda(m)*RI) * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI));
E(n,m) = Rml_prime_RI * (4 * k(n)^2 * (-1)^(m-1) * sin(k(n)*b)) / (cos(k(n)*h) * (k(n)^2 - lambda(m)^2));
    end 
end
%%%%%%%%%%%%%%%%%%%%%%%% Now write matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = zeros(N,M);
%%%%%%%% Define F00 %%%%%%%%%%
F(1,1) = (-4 * sinh(wk*b)) / (RI * log(RE/RI) * cosh(wk*h));
%%%%%% Define F0m%%%%
for m = 2:M
Rml_star_prime_RE = lambda(m)*(besseli(l, lambda(m)*RE) * dbesselk(l, lambda(m)*RI) ...
  - besselk(l, lambda(m)*RE) * dbesseli(l, lambda(m)*RI))/((besselk(l, lambda(m)*RI) ...
  * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI)));
F(1,m) = Rml_star_prime_RE * (4 * wk^2 * (-1)^(m-1) * sinh(wk*b)) / (cosh(wk*h) * (wk^2 + lambda(m)^2));
end
%%%%% Defin Fn0%%%%%%
for n = 2:N
    F(n,1) = (-4 * sin(k(n)*b)) / (RI * log(RE/RI) * cos(k(n)*h));
end
%%%%%%% Define Fnm %%%%%%%
for n = 2:N
    for m = 2:M
        Rml_star_prime_RE = lambda(m)*(besseli(l, lambda(m)*RE) * dbesselk(l, lambda(m)*RI) ...
  - besselk(l, lambda(m)*RE) * dbesseli(l, lambda(m)*RI))/((besselk(l, lambda(m)*RI) ...
  * besseli(l, lambda(m)*RE) - besselk(l, lambda(m)*RE) * besseli(l, lambda(m)*RI)));
        F(n,m) = Rml_star_prime_RE * (4 * k(n)^2 * (-1)^(m-1) * sin(k(n)*b)) / (cos(k(n)*h) * (k(n)^2 - lambda(m)^2));
    end
end
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calling the bessel0j(l,q,option) to find the positive roots of Bessel functions
 roots = bessel0j(l, Q);        % [x_1^l, x_2^l, ..., x_Q^l] (1 x Q vector)
chi = roots / RI;              % (1 x Q vector) -- only if you need chi_q^l  
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%% Write Matrix W %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = zeros(N,Q);
%%W_{0q}
for q = 1:Q
    W(1, q) = -4*wk*(chi(q)*cosh(wk*b)-wk*coth(chi(q)*b)*sinh(wk*b)) ...
              / ((chi(q)^2-wk^2) * cosh(wk*h));
end
%%% Write W_{nq}
for n = 2:N
    for q = 1:Q
        W(n, q) = -4*k(n)*( chi(q)*cos(k(n)*b) + k(n)*coth(chi(q)*b)*sin(k(n)*b)) ...
                   / ( cos(k(n)*h)*(chi(q)^2 + k(n)^2) );
    end
end
% %%%%%%%%%%%%%%%%%%%%%%%Find E1_{q}%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for q = 1:Q
     f2 = omega^2 / g;
    E1(q) = ( chi(q)*cosh(chi(q)*h) - f2*sinh(chi(q)*h) ) ...
            / ( f2*cosh(chi(q)*d) - chi(q)*sinh(chi(q)*d) );
end
G = zeros(Q, N);        % Preallocate
for q = 1:Q
    % Anonymous function for integrand
    f1 = @(r) (besselj(l, wk*r) .* besselj(l, chi(q)*r) .* r);

    % Numerical integration over [0, RI]
    int1(q, 1) = integral(f1, 0, RI);

    % G_{q0}^l formula (assuming Z0d is constant, dbesselj is your derivative of Jl)
    G(q, 1) = (-1i * a1 / omega)*Z0d* int1(q, 1) / dbesselj(l, wk*RI);
end
%% G_{qn}
for q = 1:Q
    for n = 2:N  
        % Anonymous function for integrand
    f2 = @(r) (besseli(l, k(n)*r) .* besselj(l, chi(q)*r) .* r);
    % Numerical integration over [0, RI]
    int2(q, n) = integral(f2, 0, RI);

    % G_{qn}^l formula (assuming Znd is constant, dbesselj is your derivative of Jl)
    G(q, n) = (-1i * a1 / omega)*Znd(n) * int2(q, n) / dbesseli(l, k(n)*RI);
    end
 end
% find H_{q}
 for q = 1:Q
    % Compute the diagonal value for H at index q
    H(q) = (2 / (sinh(2 * chi(q) * (h - d))) - E1(q) + 1i * a1 * chi(q) / omega ) ...
        * RI^2*( besselj(l+1, chi(q) * RI)^2 )/(2 * dbesselj(l, chi(q) * RI));
end
H_I = diag(H);   % Diagonal matrix


T_old = zeros(2*M + 2*N + Q, 1);  % Preallocate for previous unknowns
% --- Initial guess: |v_D(r)| = 0, so psi = 0

psi = zeros(Q, 1);
for iter = 1:max_iter
%Define the right hand side vector
U = zeros(2*M + 2*N + Q, 1);  % Full vector, M+N+M+N+Q elements

for m = 1:M
    if m == 1
        % Z_0^l term (first entry in Block 1)
        U(m,1) = (besselj(0, wk*RE) * sinh(wk*b)) / (b*wk * cosh(wk*h));
    else
        % Z_m^l terms (m = 2,3,...,M)
        U(m,1) = (2 * besselj(l, wk*RE) * sinh(wk*b) * wk * (-1)^(m-1)) ...
             / (b* cosh(wk*h)*(wk^2 + lambda(m)^2));
    end
end
% Block 2: Y (size N = 3)
for n = 1:N
    if n == 1
        U(n + M, 1) = -wk*dbesselj(l, wk*RE) * (2*wk*h + sinh(2*wk*h)) /(cosh(wk*h)^2);  % Y_0^l
    else
        U(n + M, 1) = 0;  % Y_n^l
    end
end
% Block 3: X (size M)
for m = 1:M
    U(m + M + N, 1) = 0;  % X_0^l, X_m^l
end

% Block 4: W (size N)
for n = 1:N
    U(n + M + N + M, 1) = 0;  % K_0^l, K_n^l
end
% Block 5: Nonlinear part
    for q = 1:Q
        U(2*M + 2*N + q) = -1i * b1 / omega * psi(q);
    end
% 

%%%%%%%%%% DEFINE THE MATRIX S%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S = [ I_M,   -A,    O1,   O2,   O3;
     -B,    I_N,   -C,   O4,   O5;
      O1,   O2,    I_M,  -D,   O3;
     -E,    O4,   -F,    I_N,  -W;
      O10,   O11,  O10,   -G,  H_I];
% Solve for the unknowns vector T (which holds [bm; an; cm; dn; eq])
 T = pinv(S) * U;   % If S may be singular or not square
  % === (C) Extract unknowns from T for current iteration ===
    bm = T(1:M);
    an = T(M+1:M+N);
    cm = T(M+N+1:2*M+N);
    dn = T(2*M+N+1:2*M+2*N);
    eq = T(2*M+2*N+1:end);
    % === (D) Update psi(q) for next iteration ===
    for q = 1:Q
        integrand = @(r) abs(v_D(N, Q, r, Z0d, wk, RI, l, dn, Znd, k, eq, chi)) ...
                        .* v_D(N, Q, r, Z0d, wk, RI, l, dn, Znd, k, eq, chi) ...
                        .* besselj(l, chi(q)*r) .* r;
        psi(q) = integral(integrand, 0, RI, 'AbsTol', 1e-8, 'RelTol', 1e-8);
     end
    if iter > 1
        diff_T = max(abs(T - T_old));
            fprintf('Iteration %d: max(|T-T_old|) = %.3e\n', iter, diff_T);
           if diff_T < tol
         converged = true;
                fprintf('  Solution converged after %d inner iterations.\n', iter);
             

            break;
            end
        end
        T_old = T;
end
 
    % --- Store results for the current 'i' (X_kRE) case ---
    % Use 'i' as the index for storing results in cell arrays.
    if converged
        bm_all{i} = bm;
        an_all{i} = an;
        cm_all{i} = cm;
        dn_all{i} = dn;
        eq_all{i} = eq;
        T_all{i}  = T;
        psi_all{i}= psi;
    else
        warning('Did not converge for X_kRE point %d. Max iterations reached.', i);
        bm_all{i} = NaN; an_all{i} = NaN; cm_all{i} = NaN; dn_all{i} = NaN; eq_all{i} = NaN;
        T_all{i}  = NaN;
        psi_all{i}= NaN;
    end
    end

Below is the helper functions which i used for my main code.

For the positive roots of bessel functions

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

For the derivative of bessel function

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

Expression for v_D(r)

function v_D_val = v_D(N, Q, r, Z0d, wk, RI, l, dn, Znd, k, eq, chi)
% Computes v_D(r) for arbitrary r, given coefficients and parameters.
% Inputs:
%   N  - number of n modes
%   Q  - number of q modes
%   r  - radial coordinate (scalar)
%   Z0, wk, RI, l, dn, Zn, k, eq, chi - see above documentation

    % First term: n=0 mode
   term1 = (dn(1) * Z0d * besselj(l, wk*r)) / (wk*dbesselj(l, wk*RI));

    % Second term: sum over n=2:N
    sum1=0;
    for nidx = 2:N % For N=3, this runs for nidx = 2 and nidx = 3
        sum1= sum1 + dn(nidx) * Znd(nidx) * besseli(l, k(nidx)*r) /(k(nidx)* dbesseli(l, k(nidx)*RI));
    end

    % Third term: sum over q=1:Q
    sum2 =0;
    for qidx = 1:Q
        sum2= sum2+ eq(qidx) * chi(qidx) * besselj(l, chi(qidx)*r) / (chi(qidx)*dbesselj(l, chi(qidx)*RI));
    end
 % Sum all terms to get the final v_D(r) value
    v_D_val = term1 + sum1 + sum2;
end

[Torsten]

The first thing that I'd do is to output the matrices A, B, C, D E, F, W, G and H_I and inspect where the entries in the order of 1e13 - 1e29 stem from in your computations.