As indicated in the: "wave propagation in structures: an FFT-based spectral analysis methodology", by James F. Doyle:
The spectral solution of 1D- elastic wave equation is as follow:
where 
is the wavenumber. C and D are the undetermined amplitudes at each frequency. Let the end of the bar at x = 0 be subjected to a force history F(t), that is, 
. E and A are the Elastic moduli and cross-sectional area respectively. The final solution is the inverse Fourier transform of the following expression:
is the Fourier transform of the applied force F(t).The numerical example for the above problem is provided as follow:
Rod:
diameter=1 inch
density=0.00247 lb/ci
E=10.6e6 lb/si
Pulse, F(t):
0.000000 0
0.001000 0
0.001100 1000
0.001300 0
0.001500 0
(sec) (N)
I wrote the following code in MatLab:
clear all
close all
clc
d=1.0;
A=pi/4*d^2;
rho=0.000247;
E=10.6e6;
n=2^15;
dt=5e-6;
fs=1/dt;
time_fcn = (0:n-1)/fs;
frequency = (0:n-1)*(fs/n);
omega=2*pi*frequency;
F=zeros(1,numel(time_fcn));
nn=find(time_fcn>=0.0011 & time_fcn<=0.0013);
F(nn)=-5e6*(time_fcn(nn)-0.0013);
plot(time_fcn,F)
Fn=fft(F);
plot(omega,Fn)
k=omega*sqrt(0.000247/10.6e6);
A=-Fn(2:numel(omega))./(1i*k(2:numel(omega))*E*A);
x=0;
G(2:numel(omega))= A.*exp(-1i*k(2:numel(omega))*x);
G(1)=simpsons(F,0,max(time_fcn),numel(time_fcn));
U= ifft(G);
plot(time_fcn*1000,U)
The result must be as follow:
However, I cannot get the same result as indicated in the abovementioned book. Can anybody tell me where is my mistake?
Thank you all,
Regards.
