how to solve a system of 12 equations using mathnet

Question

I need to solve system of 12 equations in C#. in the system phi2, phi3, wq2, wq3, eq2,eq3 are unknown variables. L1, L2, L3 are entered from the keyboard. phi1 changes from 0 to 2Pi

L1*cos(phi1)+L2*cos(phi2)-X0+L3*cos(phi3);
 -L1*sin(phi1)-L2*sin(phi2)+L3*sin(phi3)-Y0;
 -L1*sin(phi1)-L2*sin(phi2)*wq2-L3*sin(phi3)*wq3;
 -L1*cos(phi1)-L2*cos(phi2)*wq2+L3*cos(phi3)*wq3;
 -L1*cos(phi1)-(L2*cos(phi2)*wq2*wq2+L2*sin(phi3)*eq2)-(L3*cos(phi3)*wq3*wq3+L3*sin(phi3)*eq3);
 =L1*sin(phi1)-(-L2*sin(phi2)*wq2*wq2+L2*cos(phi2)*eq2)+(-L3*sin(phi3)*wq3*wq3+L3*cos(phi3)*eq3);

I have a solution for this system in Mathcad and Matlab. I tried to translate this using gpt, but the program gave incorrect valuesmathcad

Matlab:

function kusrach
global L1 L2 L3 phi1 Y0 X0
po=7800; d1=0.2; d2=0.2; d3=0.2;
X0=11.1; 
Y0=4.0;
L1=2.0; % кривогип
L2=7.225; % шатун
L3=8.05; % коромысло
phi2=45*pi/180; phi3=105*pi/180; wq2=0; wq3=0; aq2=0; aq3=0;
Res=[];
for phi1=0:pi/180:2*pi
    [r, ~, ~]=fsolve(@Ugol, [phi2 phi3 wq2 wq3 aq2 aq3]);
    Res=[Res; phi1, r];
end
Res(:,1)=Res(:,1)*180/pi;
phi1=Res(:,1);
phi2=Res(:,2); phi3=Res(:,3); wq2=Res(:,4); wq3=Res(:,5); aq2=Res(:,6); aq3=Res(:,7);
clc;
figure(1); plot(phi1, Res(:,2)*180/pi,'k--', phi1, Res(:,3)*180/pi,'b'); grid on;
legend('fi2', 'fi3'); xlabel('fi1, град'); ylabel('fi2, fi3, град');
figure(2); plot(phi1, Res(:,4),'k--', phi1, Res(:,5),'b'); grid on;
legend('wq2', 'wq3'); xlabel('fi1, град'); ylabel('wq2, wq3, рад/сек');
figure(3); plot(phi1, Res(:,6),'k--', phi1, Res(:,7),'b'); grid on;
legend('Eq2', 'Eq3'); xlabel('fi1, град'); ylabel('Eq2, Eq3, рад/сек^2');

function f=Ugol(x)
    %x(1)=fi2, x(2)=fi3 3-wq2 4-wq3 5-eq2 6-eq3 
   f(1)=L1*cos(phi1)+L2*cos(x(1))-X0+L3*cos(x(2));
   f(2)=-L1*sin(phi1)-L2*sin(x(1))+L3*sin(x(2))-Y0;
   f(3)=-L1*sin(phi1)-L2*sin(x(1))*x(3)-L3*sin(x(2))*x(4);
   f(4)=-L1*cos(phi1)-L2*cos(x(1))*x(3)+L3*cos(x(2))*x(4);
   f(5)=-L1*cos(phi1)-(L2*cos(x(1))*x(3)*x(3)+L2*sin(x(1))*x(5))-(L3*cos(x(2))*x(4)*x(4)+L3*sin(x(2))*x(6));
   f(6)=L1*sin(phi1)-(-L2*sin(x(1))*x(3)*x(3)+L2*cos(x(1))*x(5))+(-L3*sin(x(2))*x(4)*x(4)+L3*cos(x(2))*x(6));
   
phi1=input('phi1=');
phi1=phi1+1;
phi2=phi2(phi1)*180/pi; phi3=phi3(phi1)*180/pi; wq2=wq2(phi1); wq3=wq3(phi1); aq2=aq2(phi1); aq3=aq3(phi1);;
X=['phi2=',num2str(phi2),' phi3=',num2str(phi3), ' wq2=',num2str(wq2), ' wq3=',num2str(wq3),' eq2=',num2str(aq2), ' eq3=',num2str(aq3)];
disp(X);
end

code from gpt on c#

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;

namespace kursachshavnev
{
    class Program
    {
        static double L1 = 20;
        static double L2 = 85;
        static double L3 = 70;
        static double X0 = 111;
        static double Y0 = 40;
        static double phi2 = 0;
        static double phi3 = 1.411;
        static double wq2 = 0;
        static double wq3 = 0;
        static double eq2 = 0;
        static double eq3 = 0;

        static List<double[]> Res = new List<double[]>();

        // Глобальная переменная phi1
        static double phi1;

        static void Main(string[] args)
        {
            for (phi1 = 0; phi1 <= 2 * Math.PI; phi1 += Math.PI / 180)
            {
                double[] r = fsolve(Ugol, new double[] { phi2, phi3, wq2, wq3, eq2, eq3 });
                Res.Add(new double[] { phi1, r[0], r[1], r[2], r[3], r[4], r[5] });
            }

            Res.ForEach(r => Console.WriteLine($"phi1 = {r[0] * 180 / Math.PI}, phi2 = {r[1]}, phi3 = {r[2]}, wq2 = {r[3]}, wq3 = {r[4]}, eq2 = {r[5]}, eq3 = {r[6]}"));

            Console.ReadKey();
        }

        static double[] Ugol(double[] x)
        {
            double[] f = new double[6];
            // x0=phi2 x1=phi3 x2-wq2 x3-wq3 x4-eq2 x5-eq3
            f[0] = 20 * Math.Cos(phi1) + 85 * Math.Cos(x[0]) + 70 * Math.Cos(x[1]) - 111;
            f[1] = -20 * Math.Sin(phi1) - 85 * Math.Sin(x[0]) + 70 * Math.Sin(x[1]) - 40;
            f[2] = -20 * Math.Sin(phi1) - 85 * Math.Sin(x[0]) * x[2] - 70 * Math.Sin(x[1]) * x[3];
            f[3] = -20 * Math.Cos(phi1) - 85 * Math.Cos(x[0]) * x[2] + 70 * Math.Cos(x[1]) * x[3];
            f[4] = -20 * Math.Cos(phi1) - (85 * Math.Cos(x[0]) * x[2] * x[2] + 85 * Math.Sin(x[0]) * x[4]) - (70 * Math.Cos(x[1]) * x[3] * x[3] + 70 * Math.Sin(x[1]) * x[5]);
            f[5] = 20 * Math.Sin(phi1) - (-85 * Math.Sin(x[0]) * x[2] * x[2] + 85 * Math.Cos(x[0]) * x[4]) + (-70 * Math.Sin(x[1]) * x[3] * x[3] + 70 * Math.Cos(x[1]) * x[5]);

            return f;
        }

        static double[] fsolve(Func<double[], double[]> f, double[] x0)
        {
            double[] x = x0;
            double[] dx = new double[x.Length];
            double[] f0 = f(x);
            DenseMatrix J = DenseMatrix.Create(x.Length, x.Length, x.Length);

            for (int i = 0; i < x.Length; i++)
            {
                double h = 1e-3;
                for (int j = 0; j < x.Length; j++)
                {
                    x[j] += h;
                    double[] f1 = f(x);
                    J[i, j] = (f1[i] - f0[i]) / h;
                    x[j] -= h;
                }
            }

            DenseMatrix JtJ = DenseMatrix.Create(x.Length, x.Length, x.Length);
            Vector<double> JtF = Vector<double>.Build.Dense(x.Length);

            for (int i = 0; i < x.Length; i++)
            {
                for (int j = 0; j < x.Length; j++)
                {
                    JtJ[i, j] = 0;
                    for (int k = 0; k < x.Length; k++)
                    {
                        JtJ[i, j] += J[k, i] * J[k, j];
                    }
                }

                JtF[i] = 0;
                for (int k = 0; k < x.Length; k++)
                {
                    JtF[i] += J[k, i] * f0[k];
                }
            }

            DenseMatrix JtJInv = DenseMatrix.Create(x.Length, x.Length, x.Length);

            for (int i = 0; i < x.Length; i++)
            {
                for (int j = 0; j < x.Length; j++)
                {
                    JtJInv[i, j] = JtJ[i, j];
                }
            }

            for (int i = 0; i < x.Length; i++)
            {
                for (int j = 0; j < x.Length; j++)
                {
                    if (i == j)
                    {
                        JtJInv[i, j] /= JtJ[i, i];
                    }
                    else
                    {
                        JtJInv[i, j] /= JtJ[i, j];
                        for (int k = 0; k < x.Length; k++)
                        {
                            JtJInv[i, k] -= JtJInv[i, j] * JtJ[j, k];
                        }
                    }
                }
            }

            for (int i = 0; i < x.Length; i++)
            {
                for (int j = 0; j < x.Length; j++)
                {
                    dx[i] += JtJInv[i, j] * JtF[j];
                }
            }

            for (int i = 0; i < x.Length; i++)
            {
                x[i] -= dx[i];
            }

            return x;
        }
    }
}

results for phi1=0° in mathcad results mathcad

results for phi1=0° c#