K=5
M=0.1
C=0.2
 %位移 速度 加速度初值設定
d1=0;
v1=0;
a1=0;

%外荷載向量
f=zeros(1,tend/dt+1);
dt=0.002;
tend=0.4;
omg=16;
P=500;
t0=0;
f(1,1)=0;

for i=2:(tend/dt+1);
    t(1,i)=t0+dt*(i-1);
    f(1,i)=P*sin(omg*t(1,i));
end
f

f(1,i)=f1(1,i-1);
%調用NewMark-beta法求解位移 速度 加速度
[d,v,a]=Newmark(K,M,C,f,d1,v1,a1,dt,tend);

%圖形顯示
figure(1)
plot(t,d(1,:))
figure(2)
plot(t,v(1,:))
figure(3)
plot(t,a(1,:)) %加速度顯示



nodedd=d(1,:)';
nodevv=v(1,:)';
nodeaa=a(1,:)';