In order to calculate the transfer function one can assume that $f(t) = e^{s\,t}$, with $s$ representing the Laplace variable. Using separation of variables $u(t,x) = T(t)\,X(x)$ in $(1)$ yields
$$
\frac{\ddot{T}(t)}{T(t)} = \frac{\dot{k}(x)\,\dot{X}(x) + k(x)\,\ddot{X}(x)}{m(x)\,X(x)} = c. \tag{4}\label{sov_eq1}
$$
The separation of variables together with the first boundary condition $(3a)$ can be satisfied using $T(t) = e^{s\,t}$ and $X(-L) = 1$. This also implies $c = s^2$, which can be used to write the right hand side of $\eqref{sov_eq1}$ as
$$
\ddot{X}(x) = \frac{s^2\,m(x)\,X(x) - \dot{k}(x)\,\dot{X}(x)}{k(x)}. \tag{5}\label{ode_x}
$$
By using the separation of variables with $T(t) = e^{s\,t}$, the second boundary condition $(3b)$ can also be written as
$$
s\,X(L) = -\sqrt{\frac{k(L)}{m(L)}} \dot{X}(L). \tag{6}\label{ode_x_bc}
$$
The second order ordinary differential equation from $\eqref{ode_x}$ together with the constraints $X(-L) = 1$ and $\eqref{ode_x_bc}$ is equivalent to the following linear "time" varying autonomous system
$$
\left\{
\begin{align}
\dot{z}(x) &=
\begin{bmatrix}
0 & 1 \\
\frac{s^2\,m(x)}{k(x)} & \frac{-\dot{k}(x)}{k(x)}
\end{bmatrix} z(x), \\
\text{s.t.} &\ \begin{bmatrix}1 & 0\end{bmatrix} z(-L) = 1, \quad \begin{bmatrix}s & \sqrt{\frac{k(L)}{m(L)}}\end{bmatrix} z(L) = 0,
\end{align}\right. \tag{7}\label{LTV}
$$
with $z(x) = \begin{bmatrix}X(x) & \dot{X}(x)\end{bmatrix}^\top$. It can be noted that $z(L)$ is indirectly also a function of $s$ and the transfer function can be calculated with
$$
G(s) = \frac{\mathcal{L}\{u(t,L)\}(s)}{\mathcal{L}\{f(t)\}(s)} = \{\begin{bmatrix}1 & 0\end{bmatrix} z(L)\}(s). \tag{8}\label{tf}
$$
Linear time varying systems in general don't have a closed form analytical solution. Therefore, one has to resort to a numerical solution in order to solve $\eqref{LTV}$. However, that system is linear, which means that the final state $z(L)$ is linear in the initial state $z(-L)$ according to the state transition matrix $\Phi(L,-L)$, such that $z(L) = \Phi(L,-L)\,z(-L)$. The two columns of the state transition matrix can be calculated numerically by using $\begin{bmatrix}1 & 0\end{bmatrix}^\top$ and $\begin{bmatrix}0 & 1\end{bmatrix}^\top$ as initial condition for $z(-L)$ respectively. Once an approximation of $\Phi(L,-L)$ is obtained, the initial condition condition can be calculated by solving the following system of linear equations for $z(-L)$
$$
\begin{bmatrix}
\begin{bmatrix}1 & 0\end{bmatrix} \\
\begin{bmatrix}s & \sqrt{\frac{k(L)}{m(L)}}\end{bmatrix} \Phi(L,-L)
\end{bmatrix} z(-L) =
\begin{bmatrix}1 \\ 0\end{bmatrix}, \tag{9}\label{LTV_sol}
$$
and thus the value for the transfer function from $\eqref{tf}$ for a given $s$ can be obtained by solving $\eqref{LTV_sol}$ for $z(-L)$ and use $z(L) = \Phi(L,-L)\,z(-L)$.
When there is a lot of attenuation (so when $|G(s)| \ll 1$) the obtained value for $z(L)$ would be very small, in which case the obtained numerical solution might be badly conditioned. Instead, in such case it might be more accurate to solve $\eqref{LTV}$ in the reverse $x$ direction. For example by using $\hat{z}(x) = z(-x)$, such that the dynamics from $\eqref{LTV}$ can be written as
$$
\dot{\hat{z}}(x) =
\begin{bmatrix}
0 & -1 \\
\frac{-s^2\,m(-x)}{k(-x)} & \frac{\dot{k}(-x)}{k(-x)}
\end{bmatrix} \hat{z}(x), \tag{10}
$$
with $\hat{z}(-L) = z(L)$, thus the calculated state transition matrix would be $\Phi(-L,L)$ and system of linear equations from $\eqref{LTV_sol}$ would become
$$
\begin{bmatrix}
\begin{bmatrix}1 & 0\end{bmatrix} \Phi(-L,L) \\
\begin{bmatrix}s & \sqrt{\frac{k(L)}{m(L)}}\end{bmatrix}
\end{bmatrix} z(L) =
\begin{bmatrix}1 \\ 0\end{bmatrix}. \tag{11}
$$
This can be implemented in Julia using the following code:
using DifferentialEquations, ControlSystems, ControlSystemIdentification, Plots
function main(N)
Omega = exp.(range(log(1e-2), log(1e1), N))
G = [eval_TF(1im * omega) for omega = Omega]
bodeplot(FRD(Omega, G))
end
function eval_TF(s)
L = 1.0
tspan = (-L, L)
X0 = [1.0 0.0; 0.0 1.0]
p = s^2
b = [1.0 0.0]
MdivK, KxdivK = material_parameters(L)
prob1_back = ODEProblem(dx_backward!, X0[:,1], tspan, p)
prob2_back = ODEProblem(dx_backward!, X0[:,2], tspan, p)
sol1 = solve(prob1_back, Tsit5(); reltol = 1e-3, abstol = 1e-4, saveat = tspan[2])
sol2 = solve(prob2_back, Tsit5(); reltol = 1e-3, abstol = 1e-4, saveat = tspan[2])
Phi = [sol1.u[end] sol2.u[end]]
A = [b * Phi; [s*sqrt(MdivK) 1]]
if abs.(b * (A \ b'))[1] < 1
sol1 = solve(prob1_back, Tsit5(); reltol = 1e-10, abstol = 1e-11, saveat = tspan[2])
sol2 = solve(prob2_back, Tsit5(); reltol = 1e-10, abstol = 1e-11, saveat = tspan[2])
Phi = [sol1.u[end] sol2.u[end]]
A = [b * Phi; [s*sqrt(MdivK) 1]]
return (b * (A \ b'))[1]
else
prob1_forw = ODEProblem(dx_forward!, X0[:,1], tspan, p)
prob2_forw = ODEProblem(dx_forward!, X0[:,2], tspan, p)
sol1 = solve(prob1_forw, Tsit5(); reltol = 1e-10, abstol = 1e-11, saveat = tspan[2])
sol2 = solve(prob2_forw, Tsit5(); reltol = 1e-10, abstol = 1e-11, saveat = tspan[2])
Phi = [sol1.u[end] sol2.u[end]]
A = [b; [s*sqrt(MdivK) 1]*Phi]
return (b * Phi * (A \ b'))[1]
end
end
function material_parameters(x)
L = 1
B = 5 * pi / L
MdivK = 1 # Constant wave speed
KxdivK = -90 * sin(B * x)^9 # Varying impedance
return MdivK, KxdivK
end
function dx_forward!(dx,x,p,t)
MdivK, KxdivK = material_parameters(t)
dx[1] = x[2]
dx[2] = p[1] * MdivK * x[1] - KxdivK * x[2]
end
function dx_backward!(dx,x,p,t)
MdivK, KxdivK = material_parameters(-t)
dx[1] = -x[2]
dx[2] = -p[1] * MdivK * x[1] + KxdivK * x[2]
end
Which for example can generate the following Bode plot:
