How can I compare two RGB colors to see if they look the same (to human eyes). And I’d like to be able to control the amount of similarity.

Say for example 15. What is the best way to compare the colors to have more colors in the range.

I have two idea in my mind but I don’t know if they are the best options:

% (colorN(1) to colorN(3) are the RGB value of the color)
% Method 1:
if abs(double(color1(1)) – double(color2(1))) < 15 …
&& abs(double(color1(2)) – double(color2(2))) < 15 …
&& abs(double(color1(3)) – double(color2(3))) < 15
% The colors are similar
end

% Method 2
if sum(abs(double(color1(1)) – double(color2(1))), …
abs(double(color1(2)) – double(color2(2))), …
abs(double(color1(3)) – double(color2(3)))) < 15
% The colors are similar
end

% Maybe method 3 which is a combination of the two.

I’ve tried these a little bit by myself but I wasn’t happy with the results. Anyone knows a better way to compare? Does Matlab have a function for that?How can I compare two RGB colors to see if they look the same (to human eyes). And I’d like to be able to control the amount of similarity.

Say for example 15. What is the best way to compare the colors to have more colors in the range.

I have two idea in my mind but I don’t know if they are the best options:

% (colorN(1) to colorN(3) are the RGB value of the color)
% Method 1:
if abs(double(color1(1)) – double(color2(1))) < 15 …
&& abs(double(color1(2)) – double(color2(2))) < 15 …
&& abs(double(color1(3)) – double(color2(3))) < 15
% The colors are similar
end

% Method 2
if sum(abs(double(color1(1)) – double(color2(1))), …
abs(double(color1(2)) – double(color2(2))), …
abs(double(color1(3)) – double(color2(3)))) < 15
% The colors are similar
end

% Maybe method 3 which is a combination of the two.

I’ve tried these a little bit by myself but I wasn’t happy with the results. Anyone knows a better way to compare? Does Matlab have a function for that? How can I compare two RGB colors to see if they look the same (to human eyes). And I’d like to be able to control the amount of similarity.

Say for example 15. What is the best way to compare the colors to have more colors in the range.

I have two idea in my mind but I don’t know if they are the best options:

% (colorN(1) to colorN(3) are the RGB value of the color)
% Method 1:
if abs(double(color1(1)) – double(color2(1))) < 15 …
&& abs(double(color1(2)) – double(color2(2))) < 15 …
&& abs(double(color1(3)) – double(color2(3))) < 15
% The colors are similar
end

% Method 2
if sum(abs(double(color1(1)) – double(color2(1))), …
abs(double(color1(2)) – double(color2(2))), …
abs(double(color1(3)) – double(color2(3)))) < 15
% The colors are similar
end

% Maybe method 3 which is a combination of the two.

I’ve tried these a little bit by myself but I wasn’t happy with the results. Anyone knows a better way to compare? Does Matlab have a function for that? color, image processing, image, compare, rgb MATLAB Answers — New Questions

​

I’ve been trying to solve the Navier-Lame equation for thermoelasticity in a reinforced concrete beam by incorporating PDEmodel. It’s inconvinient, but I want my Young’s modulus to be a function of temperature. My procedure is solving a thermal transient problem and then taking the obtained temperature field into my c and f coefficients. The geometry is created via multicuboid, then transformed into fegeometry and fed to the thermal transient model. Then the same geometry (before conversion to fegeometry) is assigned to the general pde model (createpde(3)). The mesh is seperate between models.
The solution I get is nonsense: displacements for a 2m beam, fixed on both ends and subjected to temperature gradient of 1500K/m is around 1.2m (looks like torsion). I cannot find any mathematical errors. The attached code shows my definition of both c and f coefficients (it is just an overview). Can I use the interpolateTemperature and evaluateTemperatureGradient effectively inside my coefficients? Perhaps there is another problem with my approach or some PDE Toolbox limitations?
Thanks for any potential help on this.
Best regards
KS
%Body of the structural problem, Rt are the thermal transient reults
model = createpde(3);
model.Geometry = g;
generateMesh(model, ‘Hmax’, 0.01);

ccoeff1 = @(location,state) myfunWithAdditionalArgs1(location,state, Rt, length(tlist));
ccoeff2 = @(location,state) myfunWithAdditionalArgs2(location,state, Rt, length(tlist));
fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
fcoeff2 = @(location,state) myfunWithAdditionalArgs4(location,state, Rt, length(tlist));
gcoef = @(location,state) myfunWithAdditionalArgs5(location,state, Rt, length(tlist));

specifyCoefficients(model,’m’,0,’d’,0,’c’, ccoeff1,’a’,0,’f’,[0;0;0], ‘Cell’, 1);
specifyCoefficients(model,m=0,d=0,c=ccoeff2,a=0,f=[0;0;0], Cell=[2, 3]);

applyBoundaryCondition(model, "neumann", g=gcoef, Face=[3, 4, 5, 6]);
applyBoundaryCondition(model,"dirichlet",Face=2,u=[0,0,0]);
applyBoundaryCondition(model,"dirichlet",Face=1,u=[0,0,0]);

result = solvepde(model)

%Defining of the coefficients-just the concrete
ccoeff1 = @(location,state) myfunWithAdditionalArgs1(location,state, Rt, length(tlist));
function cmatrix1 = myfunWithAdditionalArgs1(location, state, T, id)
n1 = 45;
nr = numel(location.x);
cmatrix1 = zeros(n1, nr);
E = 30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
nu = 0.2;
mu = E./(2.*(1+nu));
lambda = E.*nu./((1-2.*nu).*(1+nu));

cmatrix1 (1,:) = 2.*mu+lambda;
cmatrix1 (3,:) = mu;
cmatrix1 (6,:) = mu;
cmatrix1 (8,:) = mu;
cmatrix1 (10,:) = lambda;
cmatrix1 (16,:) = mu;
cmatrix1 (18,:) = 2.*mu+lambda;
cmatrix1 (21,:) = mu;
cmatrix1 (24,:) = mu;
cmatrix1 (28,:) = lambda;
cmatrix1 (36,:) = mu;
cmatrix1 (38,:) = lambda;
cmatrix1 (40,:) = mu;
cmatrix1 (42,:) = mu;
cmatrix1 (45,:) = 2.*mu+lambda;
end

fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
function fmatrix1 = myfunWithAdditionalArgs3(location, state, T, id)
n1=3;
nr = numel(location.x);
fmatrix1 = zeros(n1, nr);
alpha = 1.2e-5;
nu = 0.2;
E =30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
[gradTx, gradTy, gradTz] = evaluateTemperatureGradient(T, location.x, location.y, location.z, id);
fmatrix1(1,:) = E.*alpha.*gradTx./(1-2.*nu);
fmatrix1(2,:) = E.*alpha.*gradTy./(1-2.*nu);
fmatrix1(3,:) = E.*alpha.*gradTz./(1-2.*nu);
endI’ve been trying to solve the Navier-Lame equation for thermoelasticity in a reinforced concrete beam by incorporating PDEmodel. It’s inconvinient, but I want my Young’s modulus to be a function of temperature. My procedure is solving a thermal transient problem and then taking the obtained temperature field into my c and f coefficients. The geometry is created via multicuboid, then transformed into fegeometry and fed to the thermal transient model. Then the same geometry (before conversion to fegeometry) is assigned to the general pde model (createpde(3)). The mesh is seperate between models.
The solution I get is nonsense: displacements for a 2m beam, fixed on both ends and subjected to temperature gradient of 1500K/m is around 1.2m (looks like torsion). I cannot find any mathematical errors. The attached code shows my definition of both c and f coefficients (it is just an overview). Can I use the interpolateTemperature and evaluateTemperatureGradient effectively inside my coefficients? Perhaps there is another problem with my approach or some PDE Toolbox limitations?
Thanks for any potential help on this.
Best regards
KS
%Body of the structural problem, Rt are the thermal transient reults
model = createpde(3);
model.Geometry = g;
generateMesh(model, ‘Hmax’, 0.01);

ccoeff1 = @(location,state) myfunWithAdditionalArgs1(location,state, Rt, length(tlist));
ccoeff2 = @(location,state) myfunWithAdditionalArgs2(location,state, Rt, length(tlist));
fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
fcoeff2 = @(location,state) myfunWithAdditionalArgs4(location,state, Rt, length(tlist));
gcoef = @(location,state) myfunWithAdditionalArgs5(location,state, Rt, length(tlist));

specifyCoefficients(model,’m’,0,’d’,0,’c’, ccoeff1,’a’,0,’f’,[0;0;0], ‘Cell’, 1);
specifyCoefficients(model,m=0,d=0,c=ccoeff2,a=0,f=[0;0;0], Cell=[2, 3]);

applyBoundaryCondition(model, "neumann", g=gcoef, Face=[3, 4, 5, 6]);
applyBoundaryCondition(model,"dirichlet",Face=2,u=[0,0,0]);
applyBoundaryCondition(model,"dirichlet",Face=1,u=[0,0,0]);

result = solvepde(model)

%Defining of the coefficients-just the concrete
ccoeff1 = @(location,state) myfunWithAdditionalArgs1(location,state, Rt, length(tlist));
function cmatrix1 = myfunWithAdditionalArgs1(location, state, T, id)
n1 = 45;
nr = numel(location.x);
cmatrix1 = zeros(n1, nr);
E = 30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
nu = 0.2;
mu = E./(2.*(1+nu));
lambda = E.*nu./((1-2.*nu).*(1+nu));

cmatrix1 (1,:) = 2.*mu+lambda;
cmatrix1 (3,:) = mu;
cmatrix1 (6,:) = mu;
cmatrix1 (8,:) = mu;
cmatrix1 (10,:) = lambda;
cmatrix1 (16,:) = mu;
cmatrix1 (18,:) = 2.*mu+lambda;
cmatrix1 (21,:) = mu;
cmatrix1 (24,:) = mu;
cmatrix1 (28,:) = lambda;
cmatrix1 (36,:) = mu;
cmatrix1 (38,:) = lambda;
cmatrix1 (40,:) = mu;
cmatrix1 (42,:) = mu;
cmatrix1 (45,:) = 2.*mu+lambda;
end

fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
function fmatrix1 = myfunWithAdditionalArgs3(location, state, T, id)
n1=3;
nr = numel(location.x);
fmatrix1 = zeros(n1, nr);
alpha = 1.2e-5;
nu = 0.2;
E =30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
[gradTx, gradTy, gradTz] = evaluateTemperatureGradient(T, location.x, location.y, location.z, id);
fmatrix1(1,:) = E.*alpha.*gradTx./(1-2.*nu);
fmatrix1(2,:) = E.*alpha.*gradTy./(1-2.*nu);
fmatrix1(3,:) = E.*alpha.*gradTz./(1-2.*nu);
end I’ve been trying to solve the Navier-Lame equation for thermoelasticity in a reinforced concrete beam by incorporating PDEmodel. It’s inconvinient, but I want my Young’s modulus to be a function of temperature. My procedure is solving a thermal transient problem and then taking the obtained temperature field into my c and f coefficients. The geometry is created via multicuboid, then transformed into fegeometry and fed to the thermal transient model. Then the same geometry (before conversion to fegeometry) is assigned to the general pde model (createpde(3)). The mesh is seperate between models.
The solution I get is nonsense: displacements for a 2m beam, fixed on both ends and subjected to temperature gradient of 1500K/m is around 1.2m (looks like torsion). I cannot find any mathematical errors. The attached code shows my definition of both c and f coefficients (it is just an overview). Can I use the interpolateTemperature and evaluateTemperatureGradient effectively inside my coefficients? Perhaps there is another problem with my approach or some PDE Toolbox limitations?
Thanks for any potential help on this.
Best regards
KS
%Body of the structural problem, Rt are the thermal transient reults
model = createpde(3);
model.Geometry = g;
generateMesh(model, ‘Hmax’, 0.01);

ccoeff1 = @(location,state) myfunWithAdditionalArgs1(location,state, Rt, length(tlist));
ccoeff2 = @(location,state) myfunWithAdditionalArgs2(location,state, Rt, length(tlist));
fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
fcoeff2 = @(location,state) myfunWithAdditionalArgs4(location,state, Rt, length(tlist));
gcoef = @(location,state) myfunWithAdditionalArgs5(location,state, Rt, length(tlist));

specifyCoefficients(model,’m’,0,’d’,0,’c’, ccoeff1,’a’,0,’f’,[0;0;0], ‘Cell’, 1);
specifyCoefficients(model,m=0,d=0,c=ccoeff2,a=0,f=[0;0;0], Cell=[2, 3]);

applyBoundaryCondition(model, "neumann", g=gcoef, Face=[3, 4, 5, 6]);
applyBoundaryCondition(model,"dirichlet",Face=2,u=[0,0,0]);
applyBoundaryCondition(model,"dirichlet",Face=1,u=[0,0,0]);

result = solvepde(model)

%Defining of the coefficients-just the concrete
ccoeff1 = @(location,state) myfunWithAdditionalArgs1(location,state, Rt, length(tlist));
function cmatrix1 = myfunWithAdditionalArgs1(location, state, T, id)
n1 = 45;
nr = numel(location.x);
cmatrix1 = zeros(n1, nr);
E = 30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
nu = 0.2;
mu = E./(2.*(1+nu));
lambda = E.*nu./((1-2.*nu).*(1+nu));

cmatrix1 (1,:) = 2.*mu+lambda;
cmatrix1 (3,:) = mu;
cmatrix1 (6,:) = mu;
cmatrix1 (8,:) = mu;
cmatrix1 (10,:) = lambda;
cmatrix1 (16,:) = mu;
cmatrix1 (18,:) = 2.*mu+lambda;
cmatrix1 (21,:) = mu;
cmatrix1 (24,:) = mu;
cmatrix1 (28,:) = lambda;
cmatrix1 (36,:) = mu;
cmatrix1 (38,:) = lambda;
cmatrix1 (40,:) = mu;
cmatrix1 (42,:) = mu;
cmatrix1 (45,:) = 2.*mu+lambda;
end

fcoeff1 = @(location,state) myfunWithAdditionalArgs3(location, state, Rt, length(tlist));
function fmatrix1 = myfunWithAdditionalArgs3(location, state, T, id)
n1=3;
nr = numel(location.x);
fmatrix1 = zeros(n1, nr);
alpha = 1.2e-5;
nu = 0.2;
E =30e9.*(1-0.0001*interpolateTemperature(T, location.x, location.y, location.z, id));
[gradTx, gradTy, gradTz] = evaluateTemperatureGradient(T, location.x, location.y, location.z, id);
fmatrix1(1,:) = E.*alpha.*gradTx./(1-2.*nu);
fmatrix1(2,:) = E.*alpha.*gradTy./(1-2.*nu);
fmatrix1(3,:) = E.*alpha.*gradTz./(1-2.*nu);
end pde toolbox, pdemodel MATLAB Answers — New Questions

​

How can I hide or remove the geoaxes in my plot? When I use the command:
geobasemap(‘grayland’)
everything works as expected. Indeed, both the map and the axes are visible. However, when I try to hide the axes using the following command, nothing changes:
geobasemap(‘grayland’)
gx = geoaxes;
set(gx,’Visible’,’off’)How can I hide or remove the geoaxes in my plot? When I use the command:
geobasemap(‘grayland’)
everything works as expected. Indeed, both the map and the axes are visible. However, when I try to hide the axes using the following command, nothing changes:
geobasemap(‘grayland’)
gx = geoaxes;
set(gx,’Visible’,’off’) How can I hide or remove the geoaxes in my plot? When I use the command:
geobasemap(‘grayland’)
everything works as expected. Indeed, both the map and the axes are visible. However, when I try to hide the axes using the following command, nothing changes:
geobasemap(‘grayland’)
gx = geoaxes;
set(gx,’Visible’,’off’) geobasemap, geoaxes MATLAB Answers — New Questions

​

for i=1:30
open excel INPUT file
read data fle names from Excel INPUT file
% internal loop
for j=1:6
process data
open one excel OUTPUT file
write results to excel file as a cell array along one row of excel output file
save excel file
close excel file
end
end
% have tried multiple methods, all of which failed, maybe not appropriate to MATLAB R2020b?for i=1:30
open excel INPUT file
read data fle names from Excel INPUT file
% internal loop
for j=1:6
process data
open one excel OUTPUT file
write results to excel file as a cell array along one row of excel output file
save excel file
close excel file
end
end
% have tried multiple methods, all of which failed, maybe not appropriate to MATLAB R2020b? for i=1:30
open excel INPUT file
read data fle names from Excel INPUT file
% internal loop
for j=1:6
process data
open one excel OUTPUT file
write results to excel file as a cell array along one row of excel output file
save excel file
close excel file
end
end
% have tried multiple methods, all of which failed, maybe not appropriate to MATLAB R2020b? open read write to excel file in loop MATLAB Answers — New Questions

​

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: %sn’, 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|| = %.3en’, 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|| = %.3en’, …
% 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=%2dn’, 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;
endI 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: %sn’, 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|| = %.3en’, 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|| = %.3en’, …
% 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=%2dn’, 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 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: %sn’, 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|| = %.3en’, 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|| = %.3en’, …
% 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=%2dn’, 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 bessel functions, system of equations MATLAB Answers — New Questions

​

Hi there!! I am trying to get through some of the MATLAB academy online courses, but they won’t open. The tab will open but then it just spins for hours and hours and won’t actually let me work on them. This has been happening intermittently for about a week now and I never had a problem before this. I am using an HP Ominibook 7. I already tried clearing my cache and cookies and restarting my computer, but still the issue persists. Does anybody have any ideas for what’s going on and how to fix it?Hi there!! I am trying to get through some of the MATLAB academy online courses, but they won’t open. The tab will open but then it just spins for hours and hours and won’t actually let me work on them. This has been happening intermittently for about a week now and I never had a problem before this. I am using an HP Ominibook 7. I already tried clearing my cache and cookies and restarting my computer, but still the issue persists. Does anybody have any ideas for what’s going on and how to fix it? Hi there!! I am trying to get through some of the MATLAB academy online courses, but they won’t open. The tab will open but then it just spins for hours and hours and won’t actually let me work on them. This has been happening intermittently for about a week now and I never had a problem before this. I am using an HP Ominibook 7. I already tried clearing my cache and cookies and restarting my computer, but still the issue persists. Does anybody have any ideas for what’s going on and how to fix it? matlab academy, loading course hangs MATLAB Answers — New Questions

​

Hello, how can i manage to read a video H265 in matlab on ubuntu 22.04 lts ? If you have an idea on how can I proceed I will be pleased to have your help.Hello, how can i manage to read a video H265 in matlab on ubuntu 22.04 lts ? If you have an idea on how can I proceed I will be pleased to have your help. Hello, how can i manage to read a video H265 in matlab on ubuntu 22.04 lts ? If you have an idea on how can I proceed I will be pleased to have your help. linux, h265, video MATLAB Answers — New Questions

​

Can I use the Reinforcement Learning toolbox to train a soft robot? None of the tutorials are for anything but rigid bots. Im not even sure if i can get a soft robot model into simulink as I couldnt figure it out earlierCan I use the Reinforcement Learning toolbox to train a soft robot? None of the tutorials are for anything but rigid bots. Im not even sure if i can get a soft robot model into simulink as I couldnt figure it out earlier Can I use the Reinforcement Learning toolbox to train a soft robot? None of the tutorials are for anything but rigid bots. Im not even sure if i can get a soft robot model into simulink as I couldnt figure it out earlier reinforcement learning, model, soft robot MATLAB Answers — New Questions

​

Hello,
I am having trouble using the MathWorksSimulation plugin in Unreal Engine for Unreal-Simulink cosimulation. I followed the steps detailed at the link below.
https://www.mathworks.com/help/sl3d/customize-scenes-in-unreal-engine-for-3d-simulations.html
When I open a project an Unreal Engine 5.3 and attempt to add the plugin, it prompts me to reload the scene. I hit yes to rebuild, and Unreal opens the scene is Visual Studio 2022, which is properly configured according to the Unreal Recommendations at the link below.
https://dev.epicgames.com/documentation/en-us/unreal-engine/setting-up-visual-studio-development-environment-for-cplusplus-projects-in-unreal-engine?application_version=5.3
Once I have rebuilt in Visual Studio 2022, and attempt to re-launch the project, I get this error:
"Plugin ‘MathWorksSimulation’ failed to load because module ‘MathWorksSimulation’ could not be loaded. There may be an operating system error or the module may not be properly set up."
Does anyone know how to resolve this? My computer meets the hardware and software specifications detailed here. All the files are in the correct location for Unreal to access them as far as I can verify.
Thanks in advance!
KoryHello,
I am having trouble using the MathWorksSimulation plugin in Unreal Engine for Unreal-Simulink cosimulation. I followed the steps detailed at the link below.
https://www.mathworks.com/help/sl3d/customize-scenes-in-unreal-engine-for-3d-simulations.html
When I open a project an Unreal Engine 5.3 and attempt to add the plugin, it prompts me to reload the scene. I hit yes to rebuild, and Unreal opens the scene is Visual Studio 2022, which is properly configured according to the Unreal Recommendations at the link below.
https://dev.epicgames.com/documentation/en-us/unreal-engine/setting-up-visual-studio-development-environment-for-cplusplus-projects-in-unreal-engine?application_version=5.3
Once I have rebuilt in Visual Studio 2022, and attempt to re-launch the project, I get this error:
"Plugin ‘MathWorksSimulation’ failed to load because module ‘MathWorksSimulation’ could not be loaded. There may be an operating system error or the module may not be properly set up."
Does anyone know how to resolve this? My computer meets the hardware and software specifications detailed here. All the files are in the correct location for Unreal to access them as far as I can verify.
Thanks in advance!
Kory Hello,
I am having trouble using the MathWorksSimulation plugin in Unreal Engine for Unreal-Simulink cosimulation. I followed the steps detailed at the link below.
https://www.mathworks.com/help/sl3d/customize-scenes-in-unreal-engine-for-3d-simulations.html
When I open a project an Unreal Engine 5.3 and attempt to add the plugin, it prompts me to reload the scene. I hit yes to rebuild, and Unreal opens the scene is Visual Studio 2022, which is properly configured according to the Unreal Recommendations at the link below.
https://dev.epicgames.com/documentation/en-us/unreal-engine/setting-up-visual-studio-development-environment-for-cplusplus-projects-in-unreal-engine?application_version=5.3
Once I have rebuilt in Visual Studio 2022, and attempt to re-launch the project, I get this error:
"Plugin ‘MathWorksSimulation’ failed to load because module ‘MathWorksSimulation’ could not be loaded. There may be an operating system error or the module may not be properly set up."
Does anyone know how to resolve this? My computer meets the hardware and software specifications detailed here. All the files are in the correct location for Unreal to access them as far as I can verify.
Thanks in advance!
Kory unreal engine, unreal engine plugin, mathworkssimulation MATLAB Answers — New Questions

​

Who can choose me the parameters of 13.8/220kV transformer with capacity of 150MVA or more?Who can choose me the parameters of 13.8/220kV transformer with capacity of 150MVA or more? Who can choose me the parameters of 13.8/220kV transformer with capacity of 150MVA or more? simulink, transformer MATLAB Answers — New Questions

​