%咖啡制作。
k=40;c=0;m=0.1;u0=0;u0_=0;
h=0.005;
t=0:h:0.2;%如要转载，请注明来由。与作者交流qq42644515
p=10*cos(10*t);
K=k+6*m/h^2+3*c/h;
deltaP=zeros(length(t)-1,1);
u=zeros(length(t),1);u_=zeros(length(t),1);u__=zeros(length(t),1);u__(1)=(p(1)-c*u_(1)-k*u(1))/m;
deltau=zeros(length(t)-1,1);deltau_=zeros(length(t)-1,1);deltau__=zeros(length(t)-1,1);
for i=1:length(t)-1;
    deltaP(i)=p(i+1)-p(i)+(6*m/h+3*c)*u_(i)+(3*m+c/2*h)*u__(i);
    deltau(i)=deltaP(i)/K;
    deltau_(i)=3/h*deltau(i)-3*u_(i)-u__(i)*h/2;
    deltau__(i)=6/h^2*deltau(i)-6*u_(i)/h-3*u__(i);
    u(i+1)=u(i)+deltau(i);
    u_(i+1)=u_(i)+deltau_(i);
    u__(i+1)=u__(i)+deltau__(i);
end