As stated in "Wave Propagation in Structures: an FFT-based Spectral Analysis Method" by James F. Doyle:

For the spectral solution of the 1D equation for elastic waves:

$$ widetilde {u} (x, omega) = sum Ce ^ {- i {k_ {n}} x} + De ^ {i {k_ {n}} x}, $$

from where $$ {k_ {n}} = { omega_ {n}} sqrt { rho / E} $$ is the wavenumber. C and D are the indeterminate amplitudes at each frequency. The end of the beam at x = 0 is subjected to a force curve F (t), ie$$ EA frac { partial u (x, t)} { partial x} = F (t) $$

E and A are the modulus of elasticity and the cross-sectional area, respectively. The final solution is the inverse Fourier transform of the following expression:

$$ widetilde {u} (x, omega) = – frac { widetilde {F_ {n}}} {ik_ {n} EA} e ^ {- ik_ {n} x} $$

$ widetilde {F_ {n}} $ is the Fourier transform of the applied force F (t).

The numeric example for the above problem is as follows:

```
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; %inch
A=pi/4*d^2;
rho=0.000247; %lb/inch3
E=10.6e6; %psi
%transform parameters:
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 follows:

However, I can not achieve the same result as stated in the above book. Can someone tell me where my mistake is?

Thank you all,

Greetings.