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

​

I’m trying to simulate the roll angle dynamics using the PID controller. I just would like to know if this code iw well-written or can be improved?
Any help will be appreciated.

%% Parameters
J = 0.03; % moment of inertia (kg*m^2)
Kp = 2;
Ki = 1;
Kd = 0.5;

dt = 0.01; % time step
T = 3; % total time
t = 0:dt:T;
N = length(t);

phi_ref_deg = 5;
phi_ref = deg2rad(phi_ref_deg);

phi = zeros(1, N);
omega = zeros(1, N);
u = zeros(1, N);
e = zeros(1, N);

I = 0;
prev_e = phi_ref – phi(1);

for k = 1:N-1

e(k) = phi_ref – phi(k);

I = I + e(k)*dt;
D = (e(k) – prev_e)/dt;

% PID control
u(k) = Kp*e(k) + Ki*I + Kd*D;

% Dynamics: J*phi_ddot = u
phi_ddot = u(k)/J;

omega(k+1) = omega(k) + phi_ddot*dt;
phi(k+1) = phi(k) + omega(k+1)*dt;

prev_e = e(k);
end

e(end) = phi_ref – phi(end);

phi_deg = rad2deg(phi);

peak_val = max(phi_deg);
overshoot = max(0, peak_val – phi_ref_deg);

fprintf(‘Overshoot: %.3f degn’, overshoot);

figure;
subplot(2,1,1)
plot(t, phi_deg, ‘LineWidth’, 1.4); hold on;
yline(phi_ref_deg, ‘–r’, ‘Ref’);
xlabel(‘Time (s)’);
ylabel(‘Roll angle (deg)’);
title(‘Roll Angle Response’);
grid on;

subplot(2,1,2)
plot(t, u, ‘LineWidth’, 1.4);
xlabel(‘Time (s)’);
ylabel(‘u (N*m)’);
title(‘Control Torque’);
grid on;I’m trying to simulate the roll angle dynamics using the PID controller. I just would like to know if this code iw well-written or can be improved?
Any help will be appreciated.

%% Parameters
J = 0.03; % moment of inertia (kg*m^2)
Kp = 2;
Ki = 1;
Kd = 0.5;

dt = 0.01; % time step
T = 3; % total time
t = 0:dt:T;
N = length(t);

phi_ref_deg = 5;
phi_ref = deg2rad(phi_ref_deg);

phi = zeros(1, N);
omega = zeros(1, N);
u = zeros(1, N);
e = zeros(1, N);

I = 0;
prev_e = phi_ref – phi(1);

for k = 1:N-1

e(k) = phi_ref – phi(k);

I = I + e(k)*dt;
D = (e(k) – prev_e)/dt;

% PID control
u(k) = Kp*e(k) + Ki*I + Kd*D;

% Dynamics: J*phi_ddot = u
phi_ddot = u(k)/J;

omega(k+1) = omega(k) + phi_ddot*dt;
phi(k+1) = phi(k) + omega(k+1)*dt;

prev_e = e(k);
end

e(end) = phi_ref – phi(end);

phi_deg = rad2deg(phi);

peak_val = max(phi_deg);
overshoot = max(0, peak_val – phi_ref_deg);

fprintf(‘Overshoot: %.3f degn’, overshoot);

figure;
subplot(2,1,1)
plot(t, phi_deg, ‘LineWidth’, 1.4); hold on;
yline(phi_ref_deg, ‘–r’, ‘Ref’);
xlabel(‘Time (s)’);
ylabel(‘Roll angle (deg)’);
title(‘Roll Angle Response’);
grid on;

subplot(2,1,2)
plot(t, u, ‘LineWidth’, 1.4);
xlabel(‘Time (s)’);
ylabel(‘u (N*m)’);
title(‘Control Torque’);
grid on; I’m trying to simulate the roll angle dynamics using the PID controller. I just would like to know if this code iw well-written or can be improved?
Any help will be appreciated.

%% Parameters
J = 0.03; % moment of inertia (kg*m^2)
Kp = 2;
Ki = 1;
Kd = 0.5;

dt = 0.01; % time step
T = 3; % total time
t = 0:dt:T;
N = length(t);

phi_ref_deg = 5;
phi_ref = deg2rad(phi_ref_deg);

phi = zeros(1, N);
omega = zeros(1, N);
u = zeros(1, N);
e = zeros(1, N);

I = 0;
prev_e = phi_ref – phi(1);

for k = 1:N-1

e(k) = phi_ref – phi(k);

I = I + e(k)*dt;
D = (e(k) – prev_e)/dt;

% PID control
u(k) = Kp*e(k) + Ki*I + Kd*D;

% Dynamics: J*phi_ddot = u
phi_ddot = u(k)/J;

omega(k+1) = omega(k) + phi_ddot*dt;
phi(k+1) = phi(k) + omega(k+1)*dt;

prev_e = e(k);
end

e(end) = phi_ref – phi(end);

phi_deg = rad2deg(phi);

peak_val = max(phi_deg);
overshoot = max(0, peak_val – phi_ref_deg);

fprintf(‘Overshoot: %.3f degn’, overshoot);

figure;
subplot(2,1,1)
plot(t, phi_deg, ‘LineWidth’, 1.4); hold on;
yline(phi_ref_deg, ‘–r’, ‘Ref’);
xlabel(‘Time (s)’);
ylabel(‘Roll angle (deg)’);
title(‘Roll Angle Response’);
grid on;

subplot(2,1,2)
plot(t, u, ‘LineWidth’, 1.4);
xlabel(‘Time (s)’);
ylabel(‘u (N*m)’);
title(‘Control Torque’);
grid on; code improvement 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

​

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

​

Why matlab has removed the high dimension operation of sparse?
I used to generate 3D sparse matrix using sparse. However, I find an error message when using this function says: Error using sparse, inputs must be 2-D.
Is this problem solvable?Why matlab has removed the high dimension operation of sparse?
I used to generate 3D sparse matrix using sparse. However, I find an error message when using this function says: Error using sparse, inputs must be 2-D.
Is this problem solvable? Why matlab has removed the high dimension operation of sparse?
I used to generate 3D sparse matrix using sparse. However, I find an error message when using this function says: Error using sparse, inputs must be 2-D.
Is this problem solvable? sparse matrix generation, 3d matrix MATLAB Answers — New Questions

​

I wanted to measure the THD of my circuit in MATLAB. I ran a few simulations but got absurd results. To understand, I built a test circuit to see what was happening. It observed a 0 Hz peak and an absurdly high level of THD. I’m not using any code just simscape models, so I don’t understand why this occurs. How can I fix this issue? I also attached the screenshots.I wanted to measure the THD of my circuit in MATLAB. I ran a few simulations but got absurd results. To understand, I built a test circuit to see what was happening. It observed a 0 Hz peak and an absurdly high level of THD. I’m not using any code just simscape models, so I don’t understand why this occurs. How can I fix this issue? I also attached the screenshots. I wanted to measure the THD of my circuit in MATLAB. I ran a few simulations but got absurd results. To understand, I built a test circuit to see what was happening. It observed a 0 Hz peak and an absurdly high level of THD. I’m not using any code just simscape models, so I don’t understand why this occurs. How can I fix this issue? I also attached the screenshots. fft, frequency MATLAB Answers — New Questions

​

I am interested in Architectural Acoustics, and I am wondering if you can use MATLAB with AutoCad in the design of recording studios and other spaces like concert venues, etc.?I am interested in Architectural Acoustics, and I am wondering if you can use MATLAB with AutoCad in the design of recording studios and other spaces like concert venues, etc.? I am interested in Architectural Acoustics, and I am wondering if you can use MATLAB with AutoCad in the design of recording studios and other spaces like concert venues, etc.? signal processing MATLAB Answers — New Questions

​

I have a script that is over a hundred lines long. For sake of argument, lets say lines 1 to 100 take 30 seconds to process.
I am currently working on lines 101 – 130. This 30 line block of code relies on the data that is processed from lines 1 to 100. I am new to MATLAB, so a lot of my coding from lines 101 – 130 encounters errors frequently.
The issue is I have to re-run the code from line 1 just to see if a new change I’ve made to lines 101 to 130 actually works.
I have tried using breakpoints by setting a breakpoint at line 100. But after making a change then continuing my script from line 100, I still have to run the code from line 1 if I encounter an error from line 100 to 130.
My question therefore, is can I somehow ‘cache’ the data from line 1 to 100? Is it therefore possible to run my code from line 101 to 130 multiple times over as I make constant changes through multiple errors without ever going back to line 1?I have a script that is over a hundred lines long. For sake of argument, lets say lines 1 to 100 take 30 seconds to process.
I am currently working on lines 101 – 130. This 30 line block of code relies on the data that is processed from lines 1 to 100. I am new to MATLAB, so a lot of my coding from lines 101 – 130 encounters errors frequently.
The issue is I have to re-run the code from line 1 just to see if a new change I’ve made to lines 101 to 130 actually works.
I have tried using breakpoints by setting a breakpoint at line 100. But after making a change then continuing my script from line 100, I still have to run the code from line 1 if I encounter an error from line 100 to 130.
My question therefore, is can I somehow ‘cache’ the data from line 1 to 100? Is it therefore possible to run my code from line 101 to 130 multiple times over as I make constant changes through multiple errors without ever going back to line 1? I have a script that is over a hundred lines long. For sake of argument, lets say lines 1 to 100 take 30 seconds to process.
I am currently working on lines 101 – 130. This 30 line block of code relies on the data that is processed from lines 1 to 100. I am new to MATLAB, so a lot of my coding from lines 101 – 130 encounters errors frequently.
The issue is I have to re-run the code from line 1 just to see if a new change I’ve made to lines 101 to 130 actually works.
I have tried using breakpoints by setting a breakpoint at line 100. But after making a change then continuing my script from line 100, I still have to run the code from line 1 if I encounter an error from line 100 to 130.
My question therefore, is can I somehow ‘cache’ the data from line 1 to 100? Is it therefore possible to run my code from line 101 to 130 multiple times over as I make constant changes through multiple errors without ever going back to line 1? matlab, breakpoints, cache MATLAB Answers — New Questions

​

Hello,i want to copy a plot on a uiaxes to another uiaxes and can do using the following – but I can’t get the legend to copy accross (I want to avoid using the downloadedable copyuiaxes for now)

ax1=app.UIAxes3 % My main axes for plotting
ax1.Children
L1=ax1.Legend % This works
L1s=L1.String{1} % and this

f2=figure; ax2=axes(f2);
new_children=copyobj(ax1.Children,ax2);
makeGraphBlack(app,ax2,’k’,’w’); % My own function to change appearance (background & grid)

% Copy desired axes properties
ax2.XLim=ax1.XLim; ax2.YLim=ax1.YLim;
ax2.Title.String=ax1.Title.String;
ax2.XLabel.String=ax1.XLabel.String;
ax2.YLabel.String=ax1.YLabel.String;

propsToCopy={‘FontSize’,’Color’,’GridLineStyle’};
for prop=propsToCopy
set(ax2,prop{1},get(ax1,prop{1}));
end

%Copy Legend
legend(ax2)
ax2.Legend=L1;
legend.TextColor = ‘w’;

This gives me the error message:
Unrecognized function or variable ‘legend’.

Also I notice copyobj doesn’t retain the colours of the plots, so I thought adding ‘Color’ to propsToCopy would do it. Do I have to go the route of finding the line objects of the ax1 children, and then dig down into each of their colors – or is there a quicker way?
ThanksHello,i want to copy a plot on a uiaxes to another uiaxes and can do using the following – but I can’t get the legend to copy accross (I want to avoid using the downloadedable copyuiaxes for now)

ax1=app.UIAxes3 % My main axes for plotting
ax1.Children
L1=ax1.Legend % This works
L1s=L1.String{1} % and this

f2=figure; ax2=axes(f2);
new_children=copyobj(ax1.Children,ax2);
makeGraphBlack(app,ax2,’k’,’w’); % My own function to change appearance (background & grid)

% Copy desired axes properties
ax2.XLim=ax1.XLim; ax2.YLim=ax1.YLim;
ax2.Title.String=ax1.Title.String;
ax2.XLabel.String=ax1.XLabel.String;
ax2.YLabel.String=ax1.YLabel.String;

propsToCopy={‘FontSize’,’Color’,’GridLineStyle’};
for prop=propsToCopy
set(ax2,prop{1},get(ax1,prop{1}));
end

%Copy Legend
legend(ax2)
ax2.Legend=L1;
legend.TextColor = ‘w’;

This gives me the error message:
Unrecognized function or variable ‘legend’.

Also I notice copyobj doesn’t retain the colours of the plots, so I thought adding ‘Color’ to propsToCopy would do it. Do I have to go the route of finding the line objects of the ax1 children, and then dig down into each of their colors – or is there a quicker way?
Thanks Hello,i want to copy a plot on a uiaxes to another uiaxes and can do using the following – but I can’t get the legend to copy accross (I want to avoid using the downloadedable copyuiaxes for now)

ax1=app.UIAxes3 % My main axes for plotting
ax1.Children
L1=ax1.Legend % This works
L1s=L1.String{1} % and this

f2=figure; ax2=axes(f2);
new_children=copyobj(ax1.Children,ax2);
makeGraphBlack(app,ax2,’k’,’w’); % My own function to change appearance (background & grid)

% Copy desired axes properties
ax2.XLim=ax1.XLim; ax2.YLim=ax1.YLim;
ax2.Title.String=ax1.Title.String;
ax2.XLabel.String=ax1.XLabel.String;
ax2.YLabel.String=ax1.YLabel.String;

propsToCopy={‘FontSize’,’Color’,’GridLineStyle’};
for prop=propsToCopy
set(ax2,prop{1},get(ax1,prop{1}));
end

%Copy Legend
legend(ax2)
ax2.Legend=L1;
legend.TextColor = ‘w’;

This gives me the error message:
Unrecognized function or variable ‘legend’.

Also I notice copyobj doesn’t retain the colours of the plots, so I thought adding ‘Color’ to propsToCopy would do it. Do I have to go the route of finding the line objects of the ax1 children, and then dig down into each of their colors – or is there a quicker way?
Thanks uiaxes, legend, copyobj MATLAB Answers — New Questions

​