求救~~~~matlab中fsolve解非线性方程组的问题……

polly66

是一个结构动力学中求解铁木辛可梁的问题。求解非线性代数方程组，用的fsolve函数，可是不知道为什么我算出来的数据有一部分跟正确答案一样，比如（1，15，29）得0.0457是正确的，而（9，15，28）得0.4953却不对。计算结果存在跳动点，请问各位高手 这是为什么啊！很困惑！不胜感激。程序如下：
function [x,fval,exitflag]=TIMsim1205int(I,n,s)
% inputcheck
if nargin ~= 3
    error('illegal input');
end

if ~isa(I,'numeric') || length(I) ~= 1
    error('the second input must be a scalar mumber')
end

if ~isa(n,'numeric') || length(n) ~= 1
    error('the third input must be a scalar mumber')
end

if ~isa(s,'numeric') || length(s) ~= 1
    error('the fourth input must be a scalar mumber')
end

r=5/6;
h=1;
E=1;
g=3/8*E;
p = 1;
l=s*h/12^0.5;
A=p/E*(1+E/r/g);
B=p^2/r/g/E;
C=12*p/E/(h^2);
x0=zeros(n,1);
X = zeros(n,n);
Y = zeros(n,n);
Z = zeros(n,n);
w = zeros(n,1);
v = (1:n)*pi/l;
V = v.^4*E*h^2/12/p;
ddki = zeros(1,n);
ddki(I) = 1;
dkj = zeros(n,n);
ek = zeros(1,n);

for k=1:n
    for j=1:n
        vk = v(k);
        vj = v(j);
        if k~=j,
        dkj(k,j) = -1/2*vj^2*(sin((vk-vj)*l)*vk+sin((vk-vj)*l)*vj-...
            sin((vk+vj)*l)*vk+sin((vk+vj)*l)*vj)/(vk-vj)/(vk+vj);
        else
            dkj(k,j) = -1/2*vk*(-cos(vk*l)*sin(vk*l)+vk*l);
        end
    end
end
for k=1:n
    vk = v(k);
    ek(k) = 1/2*(-cos(vk*l)*sin(vk*l)+vk*l)/vk;
end
X = diag(B*ek.*(ones(1,n)-ddki));%定义矩阵X

Y = A*dkj;%%%%%%%%%%%%%%%%定义矩阵Y
Y(1:n,I) = zeros(n,1);
for k=1:n
     Y(k,k) = A*dkj(k,k)+(2*B*V(I)-C)*ek(k);
end
Y(I,I) = B*ek(I)*V(I);

Z = A.*V(I).*dkj;%%%%%%%%%定义矩阵Z
for k=1:n
     Z(k,k) = A.*V(I).*dkj(k,k)+(B*V(I).^2+C*(V(k)-V(I))).*ek(k);
end
Z(I,I) = (A*dkj(I,I)+(2*B*V(I)-C)*ek(I))*V(I);

w = A.*V(I).*dkj(1:n,I);%%%%%%%%%%%定义常数项
w(I)=A*V(I)*dkj(I,I)+B*V(I)^2*ek(I);

options = optimset('Display','iter');
[x,fval,exitflag] = fsolve(@can,x0,options);
    function y=can(q)
        y=(q(I)^2*V(I)^2*X+Y*q(I)*V(I)+Z)*q+w;
    end
end