MATLAB: Error message in a program to solve a non-linear system of equations

'non-linear' 'system of equations' 'singular matrix'

I would like to know why I receive an error message telling me that the matrix used to solve the problem is singular to working precision and how I could fix it.
% Function to solve a non-linear system of equations
function [x]=SENL(h)
n=3;
x=(rand(n,1)); % Random column vector
tol=10^-4;     % Relative tolerance of the solution
v=zeros(n,1);  % Zeros vector intended to create a canonical vector
v(1)=1;        %with this second line.
e=(2*tol*norm(x))*v; %This is the variable we use to solve the problem.
                     % I start it with double the value of tol*norm(v) so
                     % that the 'while' loop can start
i=0;                 % Variable 'i' for the iterations
while norm(e)>tol*norm(x) && i<50 % 2 different ways of exiting the loop
    i=i+1;
    f=fun(x); %----> It evaluates the function and returns a column vector
    J=funjac(x,h,n); %----> It calculates the jacobian using incremental
                     %      quotients. By numerical methods.
    e=J\f;           % This is used to calculate e=J^(-1)*f
    x=x-e;           % It updates the solution 'x'
    norm(e)          % norm(e) gives a smaller number with every iteration.
                     % With each iteration 'x' is closer to being the
                     % solution.
end
end
function [f]=fun(x) % It evaluates a f, which is a vector-valued function
                    % with vector-valued variables
f1=x(1)^2 +4*x(2) +3*x(3);
f2=x(1)   -x(2)^2 +x(3);
f3=9*x(1) +3*x(2) -x(3)^2;
f=[f1; f2; f3];
end
function [J]=funjac(x,h,n) % funjac returns the jacobian
                           % (differential matrix of the function)
I=ones(n);
for k=1:n
    alpha(k)=max(1,abs(x(k)))*sign(x(k));
    finc=fun(x+(alpha(k)*h*I(:,k)));
    f=fun(x);
    J(:,k)=(finc-f)/(alpha(k)*h);
end
end

Best Answer

Your initial random x is drawn from rand(), which is going to be in the range (0,1) exclusive.
Inside your funjack you have
alpha(k) = max(1,abs(x(k)) * sign(x(k));
as all of your entries are in the range 0 to 1, abs(x(k)) is going to be in the range 0 to 1, and max(1,that) is going to be 1. sign() of (0,1) is going to be 1 always, because rand() never generates exactly 0. So for the first run your alpha(k) are going to be all 1.
I cannot really trace further without knowing size(h) to know what your h*I(:,k) is intended to do. I(:,k) is going to be a column vector of 1's of length n, but your n is fixed not dependent on size(h). And since your I = ones(n) is going to be 1 everywhere, I don't understand why you bother to select the k'th column. If h is p x q then since I(:,k) is going to be n x 1, for it to work then q must be the same as n, and the result would be p x 1. Multiply that by a scalar gives p x 1, and add that to the vector "x" which is n x 1, and that is only going to work if p is 1 or p is also n. That is, the dimensions would agree if h is n x n, or if h is 1 x n. Either way you get a vector n x 1 which gets passed to fun() which will promptly use exactly 3 elements of it -- which would be bad news if n were ever changed to 2.
Skipping further analysis of the function because I am tired, if you look at your output of J you will see that all three columns are the same. And that's a singular matrix.
Related Question