%咖啡制作，如要转载，请注明来由。
h=0.02;k=40;wn=20;
A=(1-cos(wn*h))/(2*k);
B=(1-cos(wn*h))/(2*k);
C=cos(wn*h);
D=sin(wn*h)/wn;
A_=sin(wn*h)/(2*k)*wn;
B_=sin(wn*h)/(2*k)*wn;
C_=-sin(wn*h)*wn;
D_=cos(wn*h);
t=0:0.02:0.4;%与笔者探讨qq:42644515
u=zeros(length(t),1);
u_=zeros(length(t),1);
p=10*cos(10*t);
for i=(1:length(t)-1);
    u(i+1)=A*p(i)+B*p(i+1)+C*u(i)+D*u_(i);
    u_(i+1)=A_*p(i)+B_*p(i+1)+C_*u(i)+D_*u_(i);
end