fminsearch and rcond with exponentials did not find expected results..i dont know if im doing something wrong or not,..please help

clc; 
    c44 = 436+9;
    e15 = -0.48;
    e11 = 7.616-11;
    u11 = 100;
    d11 = 2.6;
    p = 5700;
    %c44e = 533.34;
    e11e = 5.02e-11;

    q = 1.602e-19;
    uo = 1e20;
    %i=sqrt(-1);
    c44E = c44/e15;
    E11E = e11 / e15;
    d11E = d11 / (uo * u11);
    h = 1;
    m = uo * u11;
    m1 = p * e11 * d11;
    m2 = c44 * e11 * e15^2;
    m3 = c44 * q * uo * u11;
    m4 = uo * u11 * q;
    a = d11 * m2;

    omega_range = linspace(1, 10, 10); % Define your omega range
    results = zeros(size(omega_range));
    rcond_results = zeros(size(omega_range)); % Initialize array to store rcond results

    for jj = 1:length(omega_range)
        current_w = omega_range(jj);

        % Define the objective function for fminsearch
        objective_fn = @(k_vals) calculate_rcond(k_vals, current_w, h, c44, c44E, d11, e11, E11E, d11E, e15, m, m1, m2, m3, m4, a, p,e11e);
options=optimset('MaxIter',10000,'MaxFunEvals',1000);
        % Initial guess for [k_real, k_imag]
        initial_guess = [10, 10];

        % Use fminsearch to find the values that minimize rcond
        best_k = fminsearch(objective_fn, initial_guess,options);

        % Store the best k value for the current omega
        results(jj) = complex(best_k(1), best_k(2));

        % Calculate rcond for the best_k found
        rcond_A = calculate_rcond(best_k, current_w, h, c44, c44E, d11, e11, E11E, d11E, e15, m, m1, m2, m3, m4, a, p,e11e);

        % Store the rcond value
        rcond_results(jj) = rcond_A;
    end

    % Calculate phase velocity (omega / real(k))
   c_real = omega_range ./ real(results);

    % Plot omega vs. phase velocity
    plot(omega_range, c_real, 'LineWidth', 2);
    xlabel('\omega');
    ylabel('c_real');
    title('c_real vs. \omega');

    % Display or plot the results as needed
    disp('Results:');

Results:

    disp(results);

   1.0e+02 *

  Columns 1 through 9

   7.0653 - 0.9996i   7.0713 - 8.1675i   7.0750 - 8.1692i   7.0750 - 8.1692i   7.0750 - 8.1692i   7.0754 - 0.9971i   7.0761 - 1.0031i   7.0784 - 1.0041i   7.0784 - 1.0041i

  Column 10

   7.0766 + 1.3976i

    % Check if there are roots (rcond close to zero)
    tolerance = 1e-6; % Define a tolerance level for determining if rcond is "zero"
    is_root = rcond_results < tolerance;

    if any(is_root)
        disp('Found roots for the determinant of the matrix A.');
        disp('Corresponding omega and k values:');
        for idx = find(is_root)
            fprintf('Omega: %.2f, k: %.2f + %.2fi, rcond(A): %.2e\n', omega_range(idx), real(results(idx)), imag(results(idx)), rcond_results(idx));
        end
    else
        disp('No roots found for the determinant of the matrix A within the specified tolerance.');
    end

Found roots for the determinant of the matrix A.

Corresponding omega and k values:

Omega: 1.00, k: 706.53 + -99.96i, rcond(A): 0.00e+00
Omega: 2.00, k: 707.13 + -816.75i, rcond(A): 0.00e+00
Omega: 3.00, k: 707.50 + -816.92i, rcond(A): 0.00e+00
Omega: 4.00, k: 707.50 + -816.92i, rcond(A): 0.00e+00
Omega: 5.00, k: 707.50 + -816.92i, rcond(A): 0.00e+00
Omega: 6.00, k: 707.54 + -99.71i, rcond(A): 0.00e+00
Omega: 7.00, k: 707.61 + -100.31i, rcond(A): 0.00e+00
Omega: 8.00, k: 707.84 + -100.41i, rcond(A): 0.00e+00
Omega: 9.00, k: 707.84 + -100.41i, rcond(A): 0.00e+00
Omega: 10.00, k: 707.66 + 139.76i, rcond(A): 0.00e+00

function rcond_A = calculate_rcond(k_vals, current_w, h, c44, c44E, d11, e11, E11E, d11E, e15, m, m1, m2, m3, m4, a, p,e11e)
    k_real = k_vals(1);
    k_imag = k_vals(2);

    % Calculate intermediate values b and c
    b = m1 * current_w^2 / (k_real + k_imag*1i)^2 + current_w * m2 * 1i / (k_real + k_imag*1i)^2 - m3 / (k_real + k_imag*1i)^2;
    c = (e11 * current_w * 1i - m4) * p * current_w^2 / (k_real + k_imag*1i)^4;

    % Calculate s3 and s5 values for the current omega
    s3_val = sqrt(1+(-b-sqrt(b^2-4*a * c)) /(2*a));
    s5_val = sqrt(1+(-b+sqrt(b^2-4*a*c))/(2*a));

    % Calculate p3, p5, q3, and q5 values for the current omega
    p3_val = -e15*(s3_val^2-1)/(c44*(s3_val^2-1)+p*current_w^2/(k_real+k_imag*1i)^2);
    p5_val = -e15*(s5_val^2-1)/(c44*(s5_val^2-1) + p*current_w^2/(k_real+k_imag*1i)^2);
    q3_val = -m*(s3_val^2-1)/(d11 * (s3_val^2-1) + current_w*1i/ (k_real+k_imag*1i)^2);
    q5_val = -m*(s5_val^2-1)/(d11 * (s5_val^2-1) + current_w*1i/ (k_real+k_imag*1i)^2);

    % Define exp_kh and exp_minus_kh based on k_real and k_imag
    exp_kh = exp((k_real+1i *k_imag)*h);
    exp_minus_kh = exp(-(k_real+1i*k_imag)*h);
exp_s3=exp(s3_val*(k_real+1i*k_imag)*h);
exp_s4=exp(-s3_val*(k_real+1i*k_imag)*h);
exp_s5=exp(s5_val*(k_real+1i*k_imag)*h);
exp_s6=exp(-s5_val*(k_real+1i*k_imag)*h);
    % Define matrix A
    A = [exp_kh,exp_minus_kh,(c44E*p3_val+1)*s3_val*exp_s3,(c44E*p3_val+1)*(-s3_val)*exp_s4, ...
         (c44E*p5_val+1)*s5_val*exp_s5,(c44E*p5_val+1)*(-s5_val)*exp_s6,0,0;
         -E11E * exp_kh,E11E*exp_minus_kh, (p3_val-E11E)*s3_val * exp_s3, (p3_val - E11E) * (-s3_val) * exp_s4, ...
         (p5_val-E11E) *s5_val * exp_s5, (p5_val - E11E) * (-s5_val) * exp_s6, 0, 0;
         exp_kh,- exp_minus_kh, (1 - d11E * q3_val) * s3_val * exp_s3, (1 - d11E * q3_val) * (-s3_val) * exp_s4, ...
         (1 - d11E * q5_val) * s5_val * exp_s5, (1 - d11E * q5_val) * (-s5_val) * exp_s6, 0, 0;
         1, -1, (1 - d11E * q3_val) * s3_val, (1 - d11E * q3_val) * (-s3_val), ...
         (1 - d11E * q5_val) * s5_val, (1 - d11E * q5_val) * (-s5_val), 0, 0;
         1, -1, (c44E * p3_val + 1) * s3_val, (c44E * p3_val + 1) * (-s3_val), ...
         (c44E * p5_val + 1) * s5_val, (c44E * p5_val + 1) * (-s5_val), (-c44E * s3_val) / e15, 0;
         0, 0, p3_val, p3_val, p5_val, p5_val, -1, 0;
         -E11E, E11E, (p3_val - E11E) * s3_val, (p3_val - E11E) * (-s3_val), ...
         (p5_val - E11E) * s5_val, (p5_val - E11E) * (-s5_val), 0, e11e / e15;
         1, 1, 1, 1, 1, 1, 0, -1];

    % Compute the reciprocal condition number of the matrix A
    rcond_A = rcond(A);
end

[Torsten]

What is the aim of your optimization ? Creating a matrix that is most badly conditioned ?

[Odette]

I'm trying to find complex k which make the matrix singular..since vpasolve(det(A)) takes forever i try this approach..