for j=1:length(a)
[t,y]=ode45('kutalunwensancifang',[0:T/200:50*T],[0 0.1],[],a(j),w0);
plot(a(j),y(3000:200:end,1),'b','marker','.','markersize',5)
end
xlabel('动载荷p1/N');
ylabel('位移x/m');
for j=1:length(a)
[t,y]=ode45('kutalunwensancifang',[0:T/200:200*T],[0 0.1],[],a(j),w0);
plot(a(j),y(20000:200:end,1),'b','marker','.','markersize',5)
end