<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Duhamel Integral
Problem Data
Step1:  Natural Frequency, Damped Frequency
Step2:  Computation
Preliminaries
We chose a time step and we compute a number of constants of the integration procedure that depend on the time step
Step3:  We initialize a time variable
Step4:  We compute the load, the sines and the cosines of $\omega_D t$ and their products
Step5:  The main (and only) loop in our code, we initialize A, B and a container for saving the deflections x,
then we compute the next values of A and B, the next value of x is eventually appended to the container.
Step6:  It is necessary to plot the response.
<ASSISTANT_TASK:>
Python Code:
M = 600000
T = 0.6
z = 0.10
p0 = 400000
t0, t1, t2, t3 = 0.0, 1.0, 3.0, 6.0
Explanation: Duhamel Integral
Problem Data
End of explanation
wn = 2*np.pi/T
wd = wn*np.sqrt(1-z**2)
Explanation: Natural Frequency, Damped Frequency
End of explanation
dt = 0.05
edt = np.exp(-z*wn*dt)
fac = dt/(2*M*wd)
Explanation: Computation
Preliminaries
We chose a time step and we compute a number of constants of the integration procedure that depend on the time step
End of explanation
t = dt*np.arange(1+int(t3/dt))
Explanation: We initialize a time variable
End of explanation
p = np.where(t<=t1, p0*(t-t0)/(t1-t0), np.where(t<t2, p0*(1-(t-t1)/(t2-t1)), 0))
s = np.sin(wd*t)
c = np.cos(wd*t)
sp = s*p
cp = c*p
plt.plot(t, p/1000)
plt.xlabel('Time/s')
plt.ylabel('Force/kN')
plt.xlim((t0,t3))
plt.grid();
Explanation: We compute the load, the sines and the cosines of $\omega_D t$ and their products
End of explanation
A, B, x = 0, 0, [0]
for i, _ in enumerate(t[1:], 1):
    A = A*edt+fac*(cp[i-1]*edt+cp[i])
    B = B*edt+fac*(sp[i-1]*edt+sp[i])
    x.append(A*s[i]-B*c[i])
Explanation: The main (and only) loop in our code, we initialize A, B and a container for saving the deflections x,
then we compute the next values of A and B, the next value of x is eventually appended to the container.
End of explanation
x = np.array(x)
k = M*wn**2
Dst = p/k
plt.plot(t, x*1000)
plt.plot(t, Dst*1000)
plt.xlabel('Time/s')
plt.ylabel('Deflection/mm')
plt.xlim((t0,t3))
plt.grid()
plt.show();
Explanation: It is necessary to plot the response.
End of explanation
<END_TASK>