Does I properly translate the equations in matlab code?

lambda = 2.75; %2.18; % Example initial guess for lambda
kappa_init = 0.0371; % 0.2; %0.0269; % 0.05 % Example initial guess for kappa
theta_k_init = 0.48; %0.3748; %*0.01745; %/0.0175;% 180/10; % 0.0007778*(pi/180); % Example initial guess for theta_k in radians
R_init =  450;%  347.3092; % 328.0456; %851.6597;  %   349.4944; %362.4145; % Example initial guess for R
rout = 76.3; %76.3; %74.1783;%76.4; % Define rout

% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) ((8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2))) * ...
    (-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
    + rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
    - rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))));

B1 = @(R, kappa, lambda, theta_k, rout)((8 * R^2 * (kappa^2 - omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2) * sqrt((kappa^2 - omega_m(kappa, theta_k, lambda)^4) * (kappa^2 - omega_p(kappa, theta_k, lambda)^4))) / ...
    ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2))) * ...
    (2 / (omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2 * kappa) + rout / (kappa * omega_m(kappa, theta_k, lambda) * (omega_m(kappa, theta_k, lambda)^2 - omega_p(kappa, theta_k, lambda)^2) * besselj(1, rout * omega_m(kappa, theta_k, lambda))) ...
    - rout / (kappa * omega_p(kappa, theta_k, lambda) * (omega_m(kappa, theta_k, lambda)^2 - omega_p(kappa, theta_k, lambda)^2) * besselj(1, rout * omega_p(kappa, theta_k, lambda))));


dr = @(r_values, kappa,theta_k,R) ...
      (2 * besselj(1, r_values * omega_p(lambda, kappa, theta_k)) ./ ...
                  ((omega_m(lambda, kappa, theta_k)^2 - ...
                  omega_p(lambda, kappa, theta_k)^2) * ...
                  (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
                  omega_p(lambda, kappa, theta_k)^2))) .* ...
                 ((omega_m(lambda, kappa, theta_k)^2 * ...
                 (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
                 omega_p(lambda, kappa, theta_k) .* ...
                 (A1(R, kappa, lambda, theta_k, rout) + ...
                 (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k)^2)) ./ ...
                 (kappa^2 * omega_p(lambda, kappa, theta_k)^2 * ...
                 (rout^2 - 4 * R^2)^2)) + (B1(R, kappa, lambda, theta_k, rout) ...
                 * omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k) .* ...
                 sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
                 (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa) ...              
                 - (2 * besselj(1, r_values * omega_m(lambda, kappa, theta_k)) ./ ...
                 ((omega_m(lambda, kappa, theta_k)^2 - ...
                 omega_p(lambda, kappa, theta_k)^2) * ...
                 (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k)^2))) .* ...
                 ((omega_p(lambda, kappa, theta_k)^2 * ...
                 (kappa^2 - omega_m(lambda, kappa, theta_k)^4)) ./ ...
                  omega_m(lambda, kappa, theta_k) .* ...
                 (A1(R, kappa, lambda, theta_k, rout) + ...
                 (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k)^2)) ./ ...
                 (kappa^2 * omega_m(lambda, kappa, theta_k)^2 * ...
                 (rout^2 - 4 * R^2)^2)) + (B1(R, kappa, lambda, theta_k, rout) ...
                 * omega_m(lambda, kappa, theta_k) * ...
                 omega_p(lambda, kappa, theta_k)^2 .* ...
                 sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
                 (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa) ...  
                 - (16 * r_values * R^2 * (omega_m(lambda, kappa, theta_k)^2 + ...
                 omega_p(lambda, kappa, theta_k)^2) .* ...
                 (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k)^2)) ./ ...
                 (omega_p(lambda, kappa, theta_k)^2 * ...
                 omega_m(lambda, kappa, theta_k)^2 * ...
                 (rout^2 - 4 * R^2)^2 * (kappa^2 + ...
                 omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k)^2));
                  
dtheta = @(r_values, kappa,theta_k,R) ...
             (2 * besselj(1, r_values * omega_p(lambda, kappa, theta_k)) ./ ...
            ((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
            omega_p(lambda, kappa, theta_k)^2) * ...
            (omega_m(lambda, kappa, theta_k)^2 - ...
            omega_p(lambda, kappa, theta_k)^2))) .* ...          
            (-sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
            (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
            omega_p(lambda, kappa, theta_k) .* ... 
            (A1(R, kappa, lambda, theta_k, rout) * kappa + ...
            (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
            * omega_p(lambda, kappa, theta_k)^2)) ./ ...
            (kappa * omega_p(lambda, kappa, theta_k)^2 * ...
            (rout^2 - 4 * R^2)^2)) - B1(R, kappa, lambda, theta_k, rout) ...
            * omega_p(lambda, kappa, theta_k) * ...
            (kappa^2 - omega_m(lambda, kappa, theta_k)^4)) ...       
            + (2 * besselj(1, r_values * omega_m(lambda, kappa, theta_k)) ./ ...
            ((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
            omega_p(lambda, kappa, theta_k)^2) * ...
            (omega_m(lambda, kappa, theta_k)^2 - ...
            omega_p(lambda, kappa, theta_k)^2))) .* ...         
           (sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
           (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
            omega_m(lambda, kappa, theta_k) .* ...
           (A1(R, kappa, lambda, theta_k, rout) * kappa + ...
           (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
           omega_p(lambda, kappa, theta_k)^2)) ./ ...
           (kappa * omega_m(lambda, kappa, theta_k)^2 * ...
           (rout^2 - 4 * R^2)^2)) + B1(R, kappa, lambda, theta_k, rout) * ...
           omega_m(lambda, kappa, theta_k) * ...
           (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ...         
           + (16 * r_values * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
           * omega_p(lambda, kappa, theta_k)^2) * ...
           sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
           (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ ...
           (kappa * omega_m(lambda, kappa, theta_k)^2 * ...
           omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2 * ...
           (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
           omega_p(lambda, kappa, theta_k)^2));
           
% Plotting data 
r_values = linspace(2.16, 75.6, 100).'; % specifying the range for r_values

% Plotting
figure;

Why at the boundary both the dr and dtheta are not reaching to zero? However, the analytical soln of the equations satisfied zero boundary condition. Did I make any mistake to translate those eqns (attached here in a separate file) in maltab code? Please someone suggest me the possible error (including brackets misplaced etc).

Answer 1

[John D'Errico]

Ouch. That would be a miserable set of equations to try to code. (You have my sympathies.) I did some very quick spot checks, and I saw proper use of the dotted operators where that would be necessary, etc.

Personally, I'd probably break it down into smaller chunks, making it easier to write, read, and debug.

But a very nice tool to check your equations might be to use syms like this:

syms R kappa lambda theta_k rout

% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) ((8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2))) * ...
    (-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
    + rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
    - rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))))

A1 = function_handle with value:
    @(R,kappa,lambda,theta_k,rout)((8*R^2*(kappa^2-omega_m(lambda,kappa,theta_k)^2*omega_p(lambda,kappa,theta_k)^2))/((rout^2-4*R^2)^2*(kappa^2+omega_m(lambda,kappa,theta_k)^2*omega_p(lambda,kappa,theta_k)^2)))*(-2*(omega_m(lambda,kappa,theta_k)^2+omega_p(lambda,kappa,theta_k)^2)/(omega_m(lambda,kappa,theta_k)^2*omega_p(lambda,kappa,theta_k)^2)+rout*omega_m(lambda,kappa,theta_k)^2*(kappa^2-omega_p(lambda,kappa,theta_k)^4)/(kappa^2*omega_p(lambda,kappa,theta_k)*(omega_m(lambda,kappa,theta_k)^2-omega_p(lambda,kappa,theta_k)^2)*besselj(1,rout*omega_p(lambda,kappa,theta_k)))-rout*omega_p(lambda,kappa,theta_k)^2*(kappa^2-omega_m(lambda,kappa,theta_k)^4)/(kappa^2*omega_m(lambda,kappa,theta_k)*(omega_m(lambda,kappa,theta_k)^2-omega_p(lambda,kappa,theta_k)^2)*besselj(1,rout*omega_m(lambda,kappa,theta_k))))

omega_m(lambda,kappa,theta_k)

ans = 

omega_p(lambda,kappa,theta_k)

ans = 

A1(R,kappa,lambda,theta_k,rout)

ans = 

Now you can more easily read your code, and compare it to the equations you have. Answers formats the code to look very nice. Done at the command line, you might use pretty, or do it in a livescript to get this nice formatting.

[tuhin]

Thanks! Can you please let me know how to do it this for B1, dr and dtheta. As in matlab is difficult to get such good visible equation form. Thanks in advance.

I just want the print the eqns to compare with the attach .png file eqns (including all brackets and terms). Thus I do not want to split into small parts. Please help me with this.

[Torsten]

The representation of your equations in MATLAB and in your mathematical form differ too much that you could compare them.

[tuhin]

That is the most simplest way I could represent those eqns.

In what aspect they differ (according to you)? Did you mean the matlab representation of the mathemtical eqns are wrong here?

[Torsten]

I mean: how do you want to compare dr (e.g.) in the MATLAB representation with your mathematical description ? It doesn't help in my opinion.

Be patient and compare your MATLAB function handle and your Screenshot formulae term by term. Finally you will find possible discrepancies.

syms r_values kappa theta_k R lambda rout r_values
lambda = 2.75; %2.18; % Example initial guess for lambda
kappa_init = 0.0371; % 0.2; %0.0269; % 0.05 % Example initial guess for kappa
theta_k_init = 0.48; %0.3748; %*0.01745; %/0.0175;% 180/10; % 0.0007778*(pi/180); % Example initial guess for theta_k in radians
R_init =  450;%  347.3092; % 328.0456; %851.6597;  %   349.4944; %362.4145; % Example initial guess for R
rout = 76.3; %76.3; %74.1783;%76.4; % Define rout
% Calculate omega_m and omega_p
omega_m = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) - sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
omega_p = @(lambda, kappa, theta_k) sqrt(kappa / (2 * (lambda + 1))) * sqrt((lambda + 2) * cos(theta_k) + sqrt((lambda + 2)^2 * cos(theta_k)^2 - 4 * (lambda + 1)));
% Define A1 and B1
A1 = @(R, kappa, lambda, theta_k, rout) ((8 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2)) / ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2))) * ...
    (-2 * (omega_m(lambda, kappa, theta_k)^2 + omega_p(lambda, kappa, theta_k)^2) / (omega_m(lambda, kappa, theta_k)^2 * omega_p(lambda, kappa, theta_k)^2) ...
    + rout * omega_m(lambda, kappa, theta_k)^2 * (kappa^2 - omega_p(lambda, kappa, theta_k)^4) / (kappa^2 * omega_p(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_p(lambda, kappa, theta_k))) ...
    - rout * omega_p(lambda, kappa, theta_k)^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^4) / (kappa^2 * omega_m(lambda, kappa, theta_k) * (omega_m(lambda, kappa, theta_k)^2 - omega_p(lambda, kappa, theta_k)^2) * besselj(1, rout * omega_m(lambda, kappa, theta_k))));
B1 = @(R, kappa, lambda, theta_k, rout)((8 * R^2 * (kappa^2 - omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2) * sqrt((kappa^2 - omega_m(kappa, theta_k, lambda)^4) * (kappa^2 - omega_p(kappa, theta_k, lambda)^4))) / ...
    ((rout^2 - 4 * R^2)^2 * (kappa^2 + omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2))) * ...
    (2 / (omega_m(kappa, theta_k, lambda)^2 * omega_p(kappa, theta_k, lambda)^2 * kappa) + rout / (kappa * omega_m(kappa, theta_k, lambda) * (omega_m(kappa, theta_k, lambda)^2 - omega_p(kappa, theta_k, lambda)^2) * besselj(1, rout * omega_m(kappa, theta_k, lambda))) ...
    - rout / (kappa * omega_p(kappa, theta_k, lambda) * (omega_m(kappa, theta_k, lambda)^2 - omega_p(kappa, theta_k, lambda)^2) * besselj(1, rout * omega_p(kappa, theta_k, lambda))));
dr = @(r_values, kappa,theta_k,R) ...
      (2 * besselj(1, r_values * omega_p(lambda, kappa, theta_k)) ./ ...
                  ((omega_m(lambda, kappa, theta_k)^2 - ...
                  omega_p(lambda, kappa, theta_k)^2) * ...
                  (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
                  omega_p(lambda, kappa, theta_k)^2))) .* ...
                 ((omega_m(lambda, kappa, theta_k)^2 * ...
                 (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
                 omega_p(lambda, kappa, theta_k) .* ...
                 (A1(R, kappa, lambda, theta_k, rout) + ...
                 (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k)^2)) ./ ...
                 (kappa^2 * omega_p(lambda, kappa, theta_k)^2 * ...
                 (rout^2 - 4 * R^2)^2)) + (B1(R, kappa, lambda, theta_k, rout) ...
                 * omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k) .* ...
                 sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
                 (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa) ...              
                 - (2 * besselj(1, r_values * omega_m(lambda, kappa, theta_k)) ./ ...
                 ((omega_m(lambda, kappa, theta_k)^2 - ...
                 omega_p(lambda, kappa, theta_k)^2) * ...
                 (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k)^2))) .* ...
                 ((omega_p(lambda, kappa, theta_k)^2 * ...
                 (kappa^2 - omega_m(lambda, kappa, theta_k)^4)) ./ ...
                  omega_m(lambda, kappa, theta_k) .* ...
                 (A1(R, kappa, lambda, theta_k, rout) + ...
                 (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k)^2)) ./ ...
                 (kappa^2 * omega_m(lambda, kappa, theta_k)^2 * ...
                 (rout^2 - 4 * R^2)^2)) + (B1(R, kappa, lambda, theta_k, rout) ...
                 * omega_m(lambda, kappa, theta_k) * ...
                 omega_p(lambda, kappa, theta_k)^2 .* ...
                 sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
                 (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ kappa) ...  
                 - (16 * r_values * R^2 * (omega_m(lambda, kappa, theta_k)^2 + ...
                 omega_p(lambda, kappa, theta_k)^2) .* ...
                 (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k)^2)) ./ ...
                 (omega_p(lambda, kappa, theta_k)^2 * ...
                 omega_m(lambda, kappa, theta_k)^2 * ...
                 (rout^2 - 4 * R^2)^2 * (kappa^2 + ...
                 omega_m(lambda, kappa, theta_k)^2 * ...
                 omega_p(lambda, kappa, theta_k)^2));
dr(r_values, kappa,theta_k,R)                  

ans = 

dtheta = @(r_values, kappa,theta_k,R) ...
             (2 * besselj(1, r_values * omega_p(lambda, kappa, theta_k)) ./ ...
            ((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
            omega_p(lambda, kappa, theta_k)^2) * ...
            (omega_m(lambda, kappa, theta_k)^2 - ...
            omega_p(lambda, kappa, theta_k)^2))) .* ...          
            (-sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
            (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
            omega_p(lambda, kappa, theta_k) .* ... 
            (A1(R, kappa, lambda, theta_k, rout) * kappa + ...
            (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
            * omega_p(lambda, kappa, theta_k)^2)) ./ ...
            (kappa * omega_p(lambda, kappa, theta_k)^2 * ...
            (rout^2 - 4 * R^2)^2)) - B1(R, kappa, lambda, theta_k, rout) ...
            * omega_p(lambda, kappa, theta_k) * ...
            (kappa^2 - omega_m(lambda, kappa, theta_k)^4)) ...       
            + (2 * besselj(1, r_values * omega_m(lambda, kappa, theta_k)) ./ ...
            ((kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
            omega_p(lambda, kappa, theta_k)^2) * ...
            (omega_m(lambda, kappa, theta_k)^2 - ...
            omega_p(lambda, kappa, theta_k)^2))) .* ...         
           (sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
           (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ./ ...
            omega_m(lambda, kappa, theta_k) .* ...
           (A1(R, kappa, lambda, theta_k, rout) * kappa + ...
           (16 * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 * ...
           omega_p(lambda, kappa, theta_k)^2)) ./ ...
           (kappa * omega_m(lambda, kappa, theta_k)^2 * ...
           (rout^2 - 4 * R^2)^2)) + B1(R, kappa, lambda, theta_k, rout) * ...
           omega_m(lambda, kappa, theta_k) * ...
           (kappa^2 - omega_p(lambda, kappa, theta_k)^4)) ...         
           + (16 * r_values * R^2 * (kappa^2 - omega_m(lambda, kappa, theta_k)^2 ...
           * omega_p(lambda, kappa, theta_k)^2) * ...
           sqrt((kappa^2 - omega_m(lambda, kappa, theta_k)^4) * ...
           (kappa^2 - omega_p(lambda, kappa, theta_k)^4))) ./ ...
           (kappa * omega_m(lambda, kappa, theta_k)^2 * ...
           omega_p(lambda, kappa, theta_k)^2 * (rout^2 - 4 * R^2)^2 * ...
           (kappa^2 + omega_m(lambda, kappa, theta_k)^2 * ...
           omega_p(lambda, kappa, theta_k)^2));
           
% Plotting data 
%r_values = linspace(2.16, 75.6, 100).'; % specifying the range for r_values
% Plotting
figure;

[Steven Lord]

Rather than implementing those equations as long anonymous functions, I'd recommend implementing them as function files or local functions inside a function or script file. That way you could break some of those long expressions into smaller pieces (one per term in your equations, perhaps) and also include comments explaining how they relate to the whole and/or to the mathematical equations included in whatever text / paper / etc. you're trying to implement. [Think something like "% This is equation 1.9 on page 22"] If you do this, debugging may be easier to do (you can set a breakpoint on one line and compare what it's computing with results from a textbook or from hand calculations of that quantity.)