I want to chek my function for solving system of equaton

Hi

i write matlba function for solve a systme from a EM problem

function [Is] = current()
    [f,N,Nc,a1,a2,ra,k0,Z0,lambda ,gap_angle] = parameter();
    
    % Angular positions φ_n
    n=1:N;
    phi_n = -gap_angle+ ((2*n - 1) * pi / N);  % Compute φ_n
    
    % Initialize matrices
    A = zeros(4*N, 2*N);          % Coefficient matrix (4N equations, 2N unknowns)
    b = zeros(4*N, 1);            % Right-hand side vector
    
    % Fill the coefficient matrix A and vector b
    for m = 1:N
        for n = 1:N
            % Compute angular difference and distances R1 and R2
            phi_diff = phi_n(m) - phi_n(n);
            R1 = sqrt(ra^2 + a1^2 - 2*ra*a1*cos(phi_diff));
            R2 = sqrt(ra^2 + a2^2 - 2*ra*a2*cos(phi_diff));
            
            % Bessel functions
            H0_R1 = besselh(0, 2, k0*R1);  % H_0^(2)(k0*R1)
            H0_R2 = besselh(0, 2, k0*R2);  % H_0^(2)(k0*R2)
            H1_R1 = besselh(1, 2, k0*R1);  % H_1^(2)(k0*R1)
            H1_R2 = besselh(1, 2, k0*R2);  % H_1^(2)(k0*R2)
            
            % Equation (5): Rows 1:N
            A(m, n)=(k0*Z0/4) * H0_R1;
            
            % Equation (6): Rows N+1:2N
            A(m+N, n+N) = H0_R2;
            
            % Equation (7): Rows 2N+1:3N
            A(m+2*N, n)   = (k0*Z0/4) * H0_R1;
            A(m+2*N, n+N) = -(k0*Z0/4) * H0_R2;
            
            % Equation (8): Rows 3N+1:4N
            A(m+3*N, n)   = (1j*k0/4) * ((ra - a1*cos(phi_diff)) / R1) * H1_R1;
            A(m+3*N, n+N) = -(1j*k0/4) * ((ra - a2*cos(phi_diff)) / R2) * H1_R2;
        end
        
        % Right-hand side vector b
        b(m)       =E_inc_z(ra, phi_n(m));     % Incident E_z for Eq (5)
        b(m+N)     = 0;                        % Zero for Eq (6)
        b(m+2*N)   =E_inc_z(ra, phi_n(m));     % Incident E_z for Eq (7)
        b(m+3*N)   = H_inc_phi(ra, phi_n(m));   % Incident H_phi for Eq (8)
    end
    
    % Solve the system using linsolve
    x = linsolve(A, b);
    % Extract unknowns w_n^(1) and w_n^(2)
    w1 = x(1:N);
    w2 = x(N+1:2*N);
    Is = [w1; w2]
    % Display results
    %disp('w_n^(1):');
    %disp(w1);
    %disp('w_n^(2):');
    %disp(w2);

end    % added by Sam (Editor)

the system is in

the function is correct ?

thank

George