Convergence with fsolve?

Hi all,

I was trying to solve a system of equations using fsolve that I havent' been succesfull in doing so. So I have the following 4x4 matrix equation:

where is a diagonal positive matrix (so taking its invrese is no problem), and with are numerical sinusoidal functions with period of 120deg.

In attempting to solve for the values of θ I use the function fsolve, and I have set the number of function evaluations at a sufficiently large number so this doesn't cause the algorithm to stop prematurely. After that, at the end of the solver's attempt I get the following message which obviously means that there hasn't been any convergence to a root since the residuals are quite large, yet the steps taking are extremely small:

fsolve stopped because the relative norm of the current step, 6.019541e-13, is less than

max(options.StepTolerance^2,eps) = 1.000000e-12. However, the sum of squared function

values, r = 1.748432e+00, exceeds sqrt(options.FunctionTolerance) = 1.000000e-03.

What I have tried is:

1. Changing through all the solvers, but the answer did not improve (levenberg-marquardt seemed to be the fastest however).
2. Setting as a constant, but similarly there was no improvement
3. Running the solver through a loop with the final answer of the previous step acting as the intial condition for the next.
4. Changed the tolerances so that the algorithm doesn't take too small steps (however it still can't solve as r is always greater than sqrt(options.FunctionTolerance).

Below is my code with the .mat file attached that is needed to run the code.

% Inputs

    % Load Functions
    load('Struct_Rev00.mat')

    % Engine Data
    Nc = 9 ;                                    % Number of Cylinders
    
    % Solver Inputs
    N = length(J) ;                             % Degrees of Freedom
    Phi = 2 * pi / Nc * (0:2)' ;                % DOF Phase angle
    
% Solve System
options = optimoptions('fsolve','Display','iter-detailed',...
                       'Algorithm', 'levenberg-marquardt', ...
                       'MaxIterations', 1000, 'MaxFunEvals', 5000,...
                       'FunctionTolerance', 1e-3,'OptimalityTolerance', 1e-3,'StepTolerance', 1e-3) ;
th0 = zeros(N, 1) ;

[thsol, ~, ~,Jacobian] = fsolve(@(th) sysfcn(th, J, K, Tgf, Tm2f, Phi, N), th0, options) ;

%sysfcn(thsol, J, K, Tgf, Tm2f, Phi, N) 

% System of Equations Function 
function Fu = sysfcn(th, J, K, Tgf, Tm2f, Phi, N)

% Preallocate
Fu = zeros(N, 1) ;

% Equation RHS
Theta = wrapToPi(th(1:N-1) - Phi) ;                 % Angle argument
T_2 = [reshape(Tgf(Theta), N-1, 1) ; 0] ;           % T_2 function
T_1 = diag([reshape(Tm2f(Theta), N-1, 1) ; 0])      % T_1 function

% Equation
Fu(1:N) = (J - T_1) \ (-K * th(1:N) + T_2) ;

end

Thanks for your help in advance,

KMT.