一个变截面梁的有限元求解

这是我的 matlab程序：
E=6.02e11; %弹性模量
L=[0.02 0.18 0.015 0.335 0.39 0.226 0.16 0.64 0.089 0.03 0.66 0.66 0.701 0.2 0.2...
0.086 0.175 0.2 0.953 0.4 0.4 0.088 0.424 0.04 0.02 0.06 0.08 0.055 0.045 0.127..
0.07 0.1 0.049 0.135 0.135 0.0375 0.1975 0.1975 0.0725 0.135 0.135 0.005]; %单元长度
D=[0.225 0.24 0.225 0.2738 0.2905 0.31 0.31 0.31 0.31 0.31 0.3 0.3 0.3 0.308 0.3 0.3 0.3 0.3 0.3...
0.285 0.285 0.285 0.455 0.383 0.383 0.65 0.65 0.34 0.34 0.35 0.34 0.71 0.34 0.35 0.35...
0.42 0.42 0.42 0.35 0.35 0.35 0.35]; %单元外径
d=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; %单元内径
F=[0 0 0 0 3.6685e4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.8839e4 0 0 0 0 0]; %节点外力载荷
M=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]; %节点外加弯矩载荷
s_zc=[8 15 35 41]; %支承所在节点
q=7.85.*9.8e3.*pi.*(D.^2-d.^2)/4; %单元均布载荷
I=pi.*(D.^4-d.^4)./64; %单元截面惯性矩
n=size(L,2); %获取所划分的单元长度
for i=1:1:n
k(:,:,i)=BeamElementStiffness(E,I(i),L(i));
end %每个单元的刚度矩阵
K=zeros(2*n+2,2*n+2); %定义整体刚度矩阵
for i=1:1:n
K=BeamAssemble(K,k(:,:,i),i,i+1);
end %将单元刚度矩阵整合入整体刚度矩阵
K1=SparseMatrixChange(K,s_zc); %应用波阵法进行处理
WaiLi=BeamWaiLi(L,D,d,F,M);
n1=size(s_zc,2); %获取支承个数
WaiLi1=WaiLi';
for i=1:1:n1
WaiLi1(s_zc(i),：=[];
end %应用波阵法进行处理
u=K1\WaiLi1; %解方程
运行后提示：Warning: Matrix is singular to working precision. 然后出来的结果都是NaN啊
整体思路就是：将变截面的轴系分为很多段，轴系的重量当均布载荷处理为节点的应力和弯矩！最后求解的时候我已经加了边界条件（轴承处位移为0）了，怎么还是无法求解呢？