Plot showing gaps in numerical solution

I'm working on a numerical model involving the merging of four turbulent plumes, each originating from a circular area source located at the corners of a square. My goal is to compute and visualize the merged outer contour that encloses all four interacting plumes. I have attached the code below, the equation I am solving is at the bottom of this code. When I make the plot some gaps appear and is there a better way to plot it so the gaps disappear

%% Four-Plume Contour + Area Plotting with Automatic k_crit Detection
close all; clc; clf
global r0 k theta
%%Parameter Input
r0 = 0.25; %Plume Source Radius
%% STEP 1: Find k_crit by setting θ = π/4 and solving corresponding function, g(rho) = 0 for r0
theta = pi/4;
kL = 0*(1 - r0^2)^2; % Conservative lower bound
kU = 1.05*(1 - r0^2)^2; % Previous Theoretical maximum
k_testArray = linspace(kL,kU, 501);
rho_valid = NaN(size(k_testArray));  % Store valid roots for each k in sweep
k_crit = NaN;  % First k where g(rho)=0 has a real root
contactDetected = false; % Logical flag to stop at first contact
for jj = 1:length(k_testArray)
    k = k_testArray(jj);
    rho = mnhf_secant(@FourPlumes_Equation, [0.1 0.5], 1e-6, 0);
    if rho > 0 && isreal(rho) && isfinite(rho)
        rho_valid(jj) = rho;
        if ~contactDetected
            k_crit = k;
            contactDetected = true;
        end
    end
end
fprintf('Detected critical k value (first contact with θ = π/4): k_crit = %.8f\n\n', k_crit);

Detected critical k value (first contact with θ = π/4): k_crit = 0.87670898

%% Parameter input: k values
Nt = 10001;  % Resolution, must be an odd number
theta_array = linspace(-pi/2, pi/2, Nt);
k_array = [0.6]

k_array = 0.6000

%k_array = [0.2 0.3 0.4 0.5 0.6 0.7 k_crit 0.8 0.9 1.2 1.5];  % test range values for k
area = zeros(1,length(k_array));
%% STEP 2: Loop over k values, plot, and compute area
figure(1); hold on; box on
for jj = 1:length(k_array)
    k = k_array(jj);
    % For merged Case
    if k >= k_crit
        rho_array1 = zeros(1, Nt);
        for ii = 1:Nt
            theta = theta_array(ii);
            rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
            if rho_array1(ii) < 0
                rho_array1(ii) = NaN;
            end
        end
        [x1, y1] = pol2cart(theta_array, rho_array1);
        plot(x1, y1, '-k'); axis square; xlim([0 2]); ylim([-1 1]);
        % Plot θ = π/4 lines for visual reference
        x_diag = linspace(0, 2, Nt); y_diag = x_diag;
        plot(x_diag, y_diag, 'r--'); plot(x_diag, -y_diag, 'r--');
        % Area estimation using trapezoidal rule
        iq = find(y1 >= y_diag);  % intersection
        if ~isempty(iq)
            area(jj) = 2 * abs(trapz(x1(min(iq):(Nt-1)/2+1), y1(min(iq):(Nt-1)/2+1)));
            area(jj) = area(jj) + x1(min(iq)-1)^2;
        end
    else
        % Before Merging Case
        rho_array1 = zeros(1, Nt);
        for ii = 1:Nt
            theta = theta_array(ii);
            rho_array1(ii) = mnhf_secant(@FourPlumes_Equation, [1.5 2], 1e-6, 0);
            if rho_array1(ii) < 0
                rho_array1(ii) = NaN;
            end
        end
        [x1, y1] = pol2cart(theta_array, rho_array1);
        plot(x1, y1, '-k'); axis square;
        % Detect breakdown region and recompute in "gap"
        ij = find(isnan(rho_array1(1:round(end/2))));
        if ~isempty(ij)
            theta_min = theta_array(max(ij));
            theta_array2 = linspace(theta_min, -theta_min, Nt);
            rho_array2 = zeros(1, Nt);
            for ii = 1:Nt
                theta = theta_array2(ii);
                rho_array2(ii) = mnhf_secant(@FourPlumes_Equation, [0.1 0.9], 1e-6, 0);
                if rho_array2(ii) < 0
                    rho_array2(ii) = NaN;
                end
            end
            [x2, y2] = pol2cart(theta_array2, rho_array2);
            plot(x2, y2, 'r.');
            ik = (Nt-1)/2 + find(isnan(x2((Nt-1)/2+1:end)));
            x = [x2((Nt-1)/2+1:min(ik)-1), x1(max(ij)+1:(Nt-1)/2+1)];
            y = [y2((Nt-1)/2+1:min(ik)-1), -y1(max(ij)+1:(Nt-1)/2+1)];
            plot(x, y, 'b--'); plot(x, -y, 'b--');
            % Area from symmetric region
            area(jj) = 2 * trapz(x, y);
        end
    end
    % Plot plume source circles (n = 4)
    pos1 = [1 - r0, -r0, 2*r0, 2*r0];   % (1, 0)
    pos2 = [-1 - r0, -r0, 2*r0, 2*r0];  % (-1, 0)
    pos3 = [-r0, 1 - r0, 2*r0, 2*r0];   % (0, 1)
    pos4 = [-r0, -1 - r0, 2*r0, 2*r0];  % (0, -1)
    rectangle('Position', pos1, 'Curvature', [1 1], 'FaceColor', 'k');
    rectangle('Position', pos2, 'Curvature', [1 1], 'FaceColor', 'k');
    rectangle('Position', pos3, 'Curvature', [1 1], 'FaceColor', 'k');
    rectangle('Position', pos4, 'Curvature', [1 1], 'FaceColor', 'k');
end

%% Function: Four-Plume Equation (Li & Flynn B3) ---
function f = FourPlumes_Equation(rho)
global r0 k theta
f = (rho.^2 - 2*rho*cos(theta) + 1 - r0^2) .* ...
    (rho.^2 - 2*rho*cos(theta + pi/2) + 1 - r0^2) .* ...
    (rho.^2 - 2*rho*cos(theta + pi) + 1 - r0^2) .* ...
    (rho.^2 - 2*rho*cos(theta + 1.5*pi) + 1 - r0^2) - k^2;
end


function r=mnhf_secant(Fun,x,tol,trace)
%MNHF_SECANT finds the root of "Fun" using secant scheme.
%
% Fun - name of the external function
% x - vector of length 2, (initial guesses)
% tol - tolerance
% trace - print intermediate results
%
% Usage  mnhf_secant(@poly1,[-0.5 2.0],1e-8,1)
%        poly1 is the name of the external function.
%        [-0.5 2.0] are the initial guesses for the root.

% Check inputs.
if nargin < 4, trace = 1;  end
if nargin < 3, tol = 1e-8; end
if (length(x) ~= 2)
 error('Please provide two initial guesses')
end

f = feval(Fun,x); % Fun is assumed to accept a vector

for i = 1:100
 x3 = x(1)-f(1)*(x(2)-x(1))/(f(2)-f(1));  % Update the guess.
 f3 = feval(Fun,x3);  % Function evaluation.
 if trace, fprintf(1,'%3i %12.5f %12.5f\n', i,x3,f3); end
 if abs(f3) < tol % Check for convergenece.
  r = x3; 
  return 
 else % Reset values for x(1), x(2), f(1) and f(2).
  x(1) = x(2); f(1) = f(2); x(2) = x3; f(2) = f3;
 end
end

r = NaN;
end

[Torsten]

You can again formulate the problem as to determine the 8 roots of a polynomial of degree 8 and choose the roots that physically make sense. Didn't you succeed with your former question where the roots of a polynomial of degree 4 had to be solved for ?

[Tashmid]

Yes I could with the previous one but with this one I am not succeding.

[dpb]

Show the attempt with roots instead; it's bound to be a lot more robust solution.

Although since are looking for only real roots and within some range, it may be better to use fzero to find the root around that location. See the section <Roots of Polynomials > for discussion, example.

[Torsten]

@Tashmid

If you are still interested: This function returns the eight roots of your FourPlumes_Equation. Interestingly, you get an even polynomial of degree eight in rho - thus all the coefficients for odd powers of rho are zero.

global r0 k theta
r0 = 0.25;
k = 0.6;
theta = pi/4;
rho = roots_function()

p = 

rho = 
  -0.7996 + 0.7595i
  -0.7996 - 0.7595i
  -0.5434 + 0.4825i
  -0.5434 - 0.4825i
   0.7996 + 0.7595i
   0.7996 - 0.7595i
   0.5434 + 0.4825i
   0.5434 - 0.4825i

rho = roots_function_result_copied()

rho = 
  -0.7996 + 0.7595i
  -0.7996 - 0.7595i
  -0.5434 + 0.4825i
  -0.5434 - 0.4825i
   0.7996 + 0.7595i
   0.7996 - 0.7595i
   0.5434 + 0.4825i
   0.5434 - 0.4825i

function rho = roots_function()
  global r0 k theta
  syms R0 K Theta a1 a2 a3 a4
  p1 = poly2sym([1, a1, 1-R0^2]);
  p2 = poly2sym([1, a2, 1-R0^2]);
  p3 = poly2sym([1, a3, 1-R0^2]);
  p4 = poly2sym([1, a4, 1-R0^2]);  
  p = expand(p1*p2*p3*p4-K^2);
  p = fliplr(coeffs(p,sym('x')));
  p = simplify(subs(p,[a1 a2 a3 a4],[-2*cos(Theta),-2*cos(Theta+sym(pi)/2),-2*cos(Theta+sym(pi)),-2*cos(Theta+3*sym(pi)/2)]))
  p = double(subs(p,[R0,Theta,K],[r0,theta,k]));
  rho = roots(p);
end

or simply copying the symbolic result from roots_function to make the code run faster:

function rho = roots_function_result_copied()
  global r0 k theta
  p = [1; 0; -4*r0^2; 0; 6*r0^4-4*r0^2-2*cos(4*theta); 0; -4*r0^2*(r0^2-1)^2; 0; -k^2+(r0^2-1)^4];
  rho = roots(p);
end

[Tashmid]

Thanks @Torsten