Find input for a control system such that system reaches a desired state $x_L$ at time $t=L$. Given is the equation in state space of a control system $$\dot{x}(t)=Ax(t)+Bu(t)$$ $x(t)$ is the state vector of length $n$, $A$ is a square matrix of $n\times n$, and $B$ is a vector of length $n$. 
System is assumed to be controllable i.e. matrix $R = [B\:AB\: A^2B ...A^{n-1}B]$ is full rank. The input is given to be of the form $u(t)=B^Te^{-A^Tt}\zeta$ such that the control variable $x$ is equal to $x_L$ at time $t=L$, starting from initial $x_0$. Inserting it to the standard solution, I wrote the equation below but don't know how to proceed forward.
$$
e^{-AL}x_L-x_0=\int_{0}^{L}{e^{-At}BB^Te^{-A^Tt}\zeta dt}
$$
Vector $\zeta$ is to be found.
 A: The vector $\zeta$ is constant and can be factored out of the integral, yielding
$$
e^{-A\,L} x_L - x_0 = W(L)\,\zeta, \tag{1}
$$
with
$$
W(L) = \int_0^L e^{-A\,t} B\,B^\top e^{-A^\top t} dt. \tag{2}
$$
It can be noted that $(2)$ is related to the controllability Gramian
$$
W_c(L) = \int_0^L e^{A\,(L-t)} B\,B^\top e^{A^\top (L-t)} dt = e^{A\,L}\,W(L)\,e^{A^\top L}, \tag{3}
$$
which should be full rank for any $L>0$ if $(A,B)$ is controllable. Since the inverse of $e^{A\,L}$ can always be found with $e^{-A\,L}$ also means that $(2)$ is full rank for any $L>0$ if $(A,B)$ is controllable. Thus the solution for your $\zeta$ from $(1)$ can always be found with
$$
\zeta = W(L)^{-1} \left(e^{-A\,L} x_L - x_0\right). \tag{4}
$$
The value for $W(L)$ from $(2)$ can be obtained using
\begin{align}
-A\,W(L) - W(L)\,A^\top &= \int_0^L \left(-A\,e^{-A\,t} B\,B^\top e^{-A^\top t} - e^{-A\,t} B\,B^\top e^{-A^\top t} A^\top\right) dt, \\
&= \int_0^L \frac{d}{dt} e^{-A\,t} B\,B^\top e^{-A^\top t} dt, \\
&= \left[e^{-A\,t} B\,B^\top e^{-A^\top t}\right]_0^L, \\
&= e^{-A\,L} B\,B^\top e^{-A^\top L} - B\,B^\top,
\end{align}
which is just Lyapunov equation in $W(L)$
$$
A\,W(L) + W(L)\,A^\top + Q = 0, \tag{5}
$$
with $Q = e^{-A\,L} B\,B^\top e^{-A^\top L} - B\,B^\top$. Such equation can for example be solved using the Kronecker product and vectorization
$$
\text{vec}(W(L)) = -\left(I \otimes A + A \otimes I\right)^{-1}\,\text{vec}(Q). \tag{6}
$$

It can be noted that Lyapunov equation (or more general Sylvester equation) solvers have trouble when $A$ has eigenvalues of the imaginary axis. In such case one could instead integrate the integral from $(3)$ directly or use a control input transformation $u(t,x) = -K\,x(t) + v(t)$, with $K$ such that $A - B\,K$ does not have eigenvalues of the imaginary axis. The values for the new control input could be found with $v(t) = B^\top  e^{-\mathcal{A}^\top t} \hat{\zeta}$, using $\mathcal{A} = A-B\,K$ instead of $A$ in each equation to solve for $\hat{\zeta}$. It can be noted that different $K$ will yield a different input $u$ as a function of time, but all should drive the state to $x_L$ at $t = L$.
