我编写一段求解非线性方程组的程序，怎么也调试不下去啦，前来请教各位
不胜感激！！！！
下面是全部代码，当然QB和QA是我挑出一个比较简单的非线性方程组来进行验证用的
应景做不下去了，说雅克比矩阵接近于奇异或变质矩阵
function snow_whzh
x_y=sym('[q2;q3]');
QB=sym('[0.4*cos(q2)+0.3*cos(q3);0.4*sin(q2)+0.3*sin(q3)]');
QA=sym('[-cos(pi/4)+0.4;-sin(pi/4)+0.1]');

%q=snow_NR(Q,A,q)
x_y=local_NR(QB,QA,x_y)

function A=snow_jocabi(B,q)
%求约束方程的雅可比矩阵
[n,NA]=size(B);
m=length(q);
for i=1:n
    for j=1:m
        A(i,j)=diff(B(i),q(j));
    end
end

function a=snow_norm(H)
%计算向量或者矩阵的无穷范数
[n,m]=size(H);
b=[0];
for i=1:n
    for j=1:m
        b(i,j)=0;
    end
end
for i=1:n
    for j=1:m
        %c=abs(H(i,j));
        b(i)=b(i)+abs(H(i,j));  
    end
end
a=b(1);
for i=2:n
    if b(i-1)<b(i)
        a=b(i);
    end
end

function q=snow_NR(Q,A,y)
%用迭代法求解非线性方程组
%x=numeric(x);
%n=length(x);
%Qx=snow_jocabi(Q,x);
kk=0;
function q=inline_snow_NR(Q,A,x,kk)
    n=length(x);
    Qx=snow_jocabi(Q,x);
    for i=1:n
        u(i)=kk;
    end
    for i=1:1.0e+5
        QQ=subs(Q,x(1),u(1));
        for j=2:n
        QQ=subs(QQ,x(j),u(j));
        end
        %QQ=numeric(QQ);
        Qxx=subs(Qx,x(1),u(1));
        for j=2:n
            Qxx=subs(Qxx,x(j),u(j));
        end
        %Qxx=numeric(Qxx);
        B=QQ-A;
        u0=Qxx\B;    %求解线形方程组
        disp(u);
        for j=1:n
            u_1(j)=u(j)+u0(j);        
        end
        u=u_1; 
        disp(u);
        B=numeric(B);        %将QQ-A转换成数值型 
        b_norm=snow_norm(B);  %计算QQ-A的无穷范数
        u0=numeric(u0);         %将u0转换成数值型
        u0_norm=snow_norm(u0);   %计算u的无穷范数
        if (b_norm<1.0e-6)&(u0_norm<1.0e-6)
            %u=nemeric(u);
            q=u;
            return;
        end
        [xx,yy]=size(QQx);
        if (rank(QQx)<yy)&(b_norm>1.0e-6)
            kk=kk+1;          %如果雅可比矩阵奇异并且Q-A不等于0时
            q=inline_snow_NR(kk);  %返回到inline_snow_NR，对u重新估值
        end
               
    end
end  % end  function in_local_NR(j)

p=inline_snow_NR(Q,A,y,kk);
q=p;
end   % end function x=local_NR1(Q,A,x)


function q=local_NR(Q,A,y)
%用N-R法求解非线性方程组
%x=numeric(x);
%n=length(x);
%Qx=snow_jocabi(Q,x);
kk=1;
p=inline_snow_NR(Q,A,y,kk);
q=p;

function q=inline_snow_NR(Q,A,x,kk)
%x=numeric(x);
    n=length(x);
    Qx=snow_jocabi(Q,x);
    for i=1:n
        u(i)=kk/10;
    end
    u=u';
    for i=1:1.0e+5
        QQ=subs(Q,x(1),u(1));
        for j=2:n
        QQ=subs(QQ,x(j),u(j));
        end
        QQ=numeric(QQ);
        Qxx=subs(Qx,x(1),u(1));
        for j=2:n
            Qxx=subs(Qxx,x(j),u(j));
        end
        Qxx=numeric(Qxx);
        disp(Qxx);
        A=numeric(A);
        B=QQ-A;
        disp(B);
        u0=Qxx\B;    %求解线形方程组   
        disp('u=');
        disp(u);
        disp('u0=');
        disp(u0);        
        u=u+u0;
        %for j=1:n
         %   u(j)=u(j)+u0(j);        
         %end
        %u=u_1; 
        disp('u=');
        disp(u);
        B=numeric(B);        %将QQ-A转换成数值型 
        b_norm=snow_norm(B);  %计算QQ-A的无穷范数
        u0=numeric(u0);         %将u0转换成数值型
        u0_norm=snow_norm(u0);   %计算u的无穷范数
        if (b_norm<1.0e-6)&(u0_norm<1.0e-6)
            %u=nemeric(u);
            q=u;
            return;
        end
        Qxx=numeric(Qxx);
        %[xx,yy]=size(Qxx);
        if (det(Qxx)<1.0e-16)&(b_norm>1.0e-6)
            kk=kk+1;          %如果雅可比矩阵奇异并且Q-A不等于0时
            q=inline_snow_NR(Q,A,x,kk);  %返回到inline_snow_NR，对u重新估值
        end
               
    end
    %end  % end  function in_local_NR(j)