matlab 纽马克 激励,用Newmark方法计算系统的动力学响应的matlab程序

请大家帮忙看看这个程序有什么问题？用Newmark方法计算系统的动力学响应，结果大的惊人。

function[Q,V,AA]=newmarkb

E=2.1e11;P=7850;D1=0.405;d1=0.375;D2=0.375;d2=0.335;D3=0.335;d3=0.285;D4=0.285;d4=0.225;D5=0.225;d5=0.150;

A=(pi*(D1^2-d1^2))/4;

I=(pi*(D1^4-d1^4))/64;

M1= Mass (P,A,I,0,0,13,0);

A=(pi*(D2^2-d2^2))/4;

I=(pi*(D2^4-d2^4))/64;

M2= Mass (P,A,I,13,0,26,0);

A=(pi*(D3^2-d3^2))/4;

I=(pi*(D3^4-d3^4))/64;

M3= Mass (P,A,I,26,0,39,0);

A=(pi*(D4^2-d4^2))/4;

I=(pi*(D4^4-d4^4))/64;

M4= Mass (P,A,I,39,0,52,0);

A=(pi*(D5^2-d5^2))/4;

I=(pi*(D5^4-d5^4))/64;

M5= Mass (P,A,I,52,0,65,0);

M=zeros(18,18);

M= MAssemble(M,M1,1,2);

M= MAssemble(M,M2,2,3);

M= MAssemble(M,M3,3,4);