Error in soln. Please help to solve the issue.

% Define initial parameters
lambda_init = 1.36;
R_init = 500;
rout = 76.4;

% Define c
c = sqrt(4 - (rout/R_init)^2);

% Define k_e and k_o ranges
kappa_e_values = linspace(0.0, 0.05, 20);
kappa_o_values = linspace(0.0, 0.05, 20);

% Initialize arrays to store the results
max_rd_dr = zeros(length(kappa_e_values), length(kappa_o_values));
min_rd_dr = zeros(length(kappa_e_values), length(kappa_o_values));
amp_rd_dr = zeros(length(kappa_e_values), length(kappa_o_values));
mod_rd_dr = zeros(length(kappa_e_values), length(kappa_o_values));

max_rd_dtheta = zeros(length(kappa_e_values), length(kappa_o_values));
min_rd_dtheta = zeros(length(kappa_e_values), length(kappa_o_values));
amp_rd_dtheta = zeros(length(kappa_e_values), length(kappa_o_values));
mod_rd_dtheta = zeros(length(kappa_e_values), length(kappa_o_values));

% Loop over kappa_e and kappa_o values
for i = 1:length(kappa_e_values)
    for j = 1:length(kappa_o_values)
        kappa_e = kappa_e_values(i);
        kappa_o = kappa_o_values(j);
        
        % Calculate kappa and theta_k
        kappa = sqrt(kappa_e^2 + kappa_o^2);
        theta_k = atan(kappa_o / kappa_e);
        
        % Initial guess for the solution
        solinit = bvpinit(linspace(0.0001, rout, 100), [0, 0, 0, 0]);
        
        % Solve the BVP
        sol = bvp4c(@(r, y) odefun(r, y, lambda_init, kappa, theta_k, c, R_init), @bcfun, solinit);
        
        % Extract solutions
        r = linspace(0.0001, rout, 100);
        y = deval(sol, r);
        dr_sol = y(1,:);
        dtheta_sol = y(3,:);
        
        % Calculate r*dr and r*dtheta
        r_dr = r .* dr_sol;
        r_dtheta = r .* dtheta_sol;
        
        % Calculate max, min, amplitude, and modulus of r*dr
        max_rd_dr(i, j) = max(r_dr);
        min_rd_dr(i, j) = min(r_dr);
        amp_rd_dr(i, j) = max_rd_dr(i, j) - min_rd_dr(i, j);
        mod_rd_dr(i, j) = max(abs(max_rd_dr(i, j)), abs(min_rd_dr(i, j)));
        
        % Calculate max, min, amplitude, and modulus of r*dtheta
        max_rd_dtheta(i, j) = max(r_dtheta);
        min_rd_dtheta(i, j) = min(r_dtheta);
        amp_rd_dtheta(i, j) = max_rd_dtheta(i, j) - min_rd_dtheta(i, j);
        mod_rd_dtheta(i, j) = max(abs(max_rd_dtheta(i, j)), abs(min_rd_dtheta(i, j)));
    end
end

% Define a custom colormap 'zet' that emphasizes peak values
zet_colormap = summer;

% Plot the results with kappa_e and kappa_o
subplot(1,2,1);
surf(kappa_e_values, kappa_o_values, amp_rd_dr');
xlabel('k_e');
ylabel('k_o');
zlabel('Amplitude of (r*d_r)');
title('Amplitude of r*d_r');
cb = colorbar;
colormap(zet_colormap);
caxis([min(amp_rd_dr(:)), max(amp_rd_dr(:))]);

subplot(1,2,2);
surf(kappa_e_values, kappa_o_values, amp_rd_dtheta');
xlabel('k_e');
ylabel('k_o');
zlabel('Amplitude of (r*d_\theta)');
title('Amplitude of r*d_\theta');
cb = colorbar;
colormap(zet_colormap);
caxis([min(amp_rd_dtheta(:)), max(amp_rd_dtheta(:))]);

% Define the function for the differential equations
function dydr = odefun(r, y, lambda_init, kappa, theta_k, c, R_init)
    dydr = zeros(4,1);
    dydr(1) = y(2);
    dydr(2) = -((lambda_init+1)*y(2)+1/r*(kappa*r^2*cos(theta_k)-(lambda_init+1))*y(1)-kappa*r*y(3)*sin(theta_k)+16*lambda_init*r^2/(c^4*R_init^2))/(r*(lambda_init+1));
    dydr(3) = y(4);   
    dydr(4) = -(y(4)+1/r*(kappa*r^2*cos(theta_k)-1)*y(3)+kappa*r*y(1)*sin(theta_k))/r;
end

% Boundary conditions
function res = bcfun(ya, yb)
    res = [ya(1); ya(3); yb(1); yb(3)];
end

[Star Strider]

Looking at the ‘Extra_Parameters’ vector (that I added just before the bvp4c call):

Extra_Parameters = [lambda_init, kappa, theta_k, c, R_init]

‘theta_k’ is NaN, and displaying ‘dydr’ revealed its elements always to be either 0 or NaN.

[Walter Roberson]

I suggest you replace

        theta_k = atan(kappa_o / kappa_e);

with

        theta_k = atan2(kappa_o, kappa_e);