lsqcurvefit Implementation for array xdata, matrix ydata, and symmetric matrix condition

I am looking for (guidance/solution/alternative solution) to use lsqcurvefit. Is it feasable to implement the following case? The following text is divided to two sections. Section I discusses the approach and section II provides the code (not working) to implement this case.

SECTION I

According to the documentation, I plan to use the following format

x = lsqcurvefit(fun, x0, xdata, ydata, lb, ub, options)

where

- fun is the function provided below

w is an array of length j

R(w) and L(w) are real symmetric matrices of size n x n that vary with w.

Rb and Lb are known matrices of the same size as R(w) and L(w). Rb and Lb are the equal to R(1) and L(1) respectively.

The series Rfk and Lfk are positive and real values.

Mk is a symmetric real matrix of size n x n.

If the diagonal elements of Lb are denoted as Lb_1, Lb_2, Lb_n, Mk can also be written as

 Mk = [k11_k*sqrt(Lb_1*Lfk), k12_k *sqrt(Lb_1*Lfk), k13_k *sqrt(Lb_1*Lfk),…. k1n_k *sqrt(Lb_1*Lfk);
k21_k*sqrt(Lb_2*Lfk), k22_k*sqrt(Lb_2*Lfk), k23_k *sqrt(Lb_2*Lfk),…. K2n_k *sqrt(Lb_2*Lfk);
…….
kn1_k *sqrt(Lb_n*Lfk), kn2_k *sqrt(Lb_n*Lfk), kn3_k *sqrt(Lb_n*Lfk),…. Knn_k *sqrt(Lb_n*Lfk);]

where Lfk are the same as the ones shown in the series, and k11_k …knn_k range from -1 to 1 and k ranges from 1 to r.

R(w) and L(w) can be related using the following equation Z = R(w) + 1i*w*L(w).

-x0 is a matrix that consists of column vectors of [(Rf1;…;Rfk), (Lf1;…;Lfk) , (k11_1; …;knn_1, k11_2; …;knn_2, k11_k; …;knn_k)]

-xdata is w

-ydata is Z

- do not know how to implement the symmetric matrix condition Mk

options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');

SECTION II

clear
clc
f = readmatrix("freq.txt");
L11 = readmatrix("L11.txt");
L12 = readmatrix("L12.txt");
L13 = readmatrix("L13.txt");
L14 = readmatrix("L14.txt");
L22 = readmatrix("L22.txt");
L23 = readmatrix("L23.txt");
L24 = readmatrix("L24.txt");
L33 = readmatrix("L33.txt");
L34 = readmatrix("L34.txt");
L44 = readmatrix("L44.txt");

R11 = readmatrix("R11.txt");
R12 = readmatrix("R12.txt");
R13 = readmatrix("R13.txt");
R14 = readmatrix("R14.txt");
R22 = readmatrix("R22.txt");
R23 = readmatrix("R23.txt");
R24 = readmatrix("R24.txt");
R33 = readmatrix("R33.txt");
R34 = readmatrix("R34.txt");
R44 = readmatrix("R44.txt");

w = 2*pi*f;

n = 4;  
L = zeros(n,n,length(f)); 
R = L;
Z = L;
m = 6;

Lf = eye(m); 
Rf = Lf;    

for i = 1:length(f)
    L(:,:,i) = [L11(i), L12(i), L13(i), L14(i);
        L12(i), L22(i), L23(i), L24(i);
        L13(i), L23(i), L33(i), L34(i);
        L14(i), L24(i), L34(i), L44(i);
        ];
    
    R(:,:,i) = [R11(i), R12(i), R13(i), R14(i);
        R12(i), R22(i), R23(i), R24(i);
        R13(i), R23(i), R33(i), R34(i);
        R14(i), R24(i), R34(i), R44(i);
        ];
    
    Z(:,:,i) = R(:,:,i) + 1j*w(i).*L(:,:,i); 
    
end

ydata = Z;

for i = 1:length(f)
xdata(:,:,i) = repelem(w(i),n,n);
end

R_f_guess = [1, 0.1, 1, 0.1]';
L_f_guess = [1, 0.1, 1, 0.1]'*1e-6;
k_guess = [ones(4,16)*0.05];
x0 = [R_f_guess, L_f_guess, k_guess];

x = lsqcurvefit(@RL,x0,xdata,ydata);


function Z = RL(x,xdata)

R_f = x(:,1)';
L_f = x(:,2)';
k = x(:,3:end)';

w = xdata(1,1,:);

sumR = 0;
sumL = 0;

r = 4;

 Lb =   1.0e-04 *[0.1063,    0.1063,    0.1063,    0.1063;
    0.1063,    0.1063,    0.1063,    0.1063;
    0.1063,    0.1063,    0.1063,    0.1063;
    0.1063,    0.1063,    0.1063,    0.1063;
];

for k = 1:r
    
    M = [k(1,r)*sqrt(Lb(1,1,1)*L_f(r)),k(2,r)*sqrt(Lb(1,1,1)*L_f(r)),k(3,r)*sqrt(Lb(1,1,1)*L_f(r)),k(4,r)*sqrt(Lb(1,1,1)*L_f(r));
    k(5,r)*sqrt(Lb(2,2,1)*L_f(r)),k(6,r)*sqrt(Lb(2,2,1)*L_f(r)),k(7,r)*sqrt(Lb(2,2,1)*L_f(r)),k(8,r)*sqrt(Lb(2,2,1)*L_f(r));
    k(9,r)*sqrt(Lb(3,3,1)*L_f(r)),k(10,r)*sqrt(Lb(3,3,1)*L_f(r)),k(11,r)*sqrt(Lb(3,3,1)*L_f(r)),k(12,r)*sqrt(Lb(3,3,1)*L_f(r));
    k(13,r)*sqrt(Lb(4,4,1)*L_f(r)),k(14,r)*sqrt(Lb(4,4,1)*L_f(r)),k(15,r)*sqrt(Lb(4,4,1)*L_f(r)),k(16,r)*sqrt(Lb(4,4,1)*L_f(r));];
    
    sumR = sumR + w.^2 * R_f(1,r) ./ (R_f(1,r).^2 + w.^2 * L_f(1,r).^2) .* M*transpose(M);
    sumL = sumL + w.^2 * L_f(1,r) ./ (R_f(1,r).^2 + w.^2 * L_f(1,r).^2) .* M*transpose(M);
end
    
    Z = Rb + sumR + 1j*w.*(Lb - sumL);  
    

end