The solution was
$$
u(x,t) = \int\limits_0^t u(x,t;s) \, ds \quad x \in \mathbb{R}^n, t \ge 0
$$
Partial differentiation for $t$ of the integral function according to Leibniz
$$
\frac{\partial}{\partial t} \int\limits_0^{b(t)} u(x, t; s) \, ds =
u(x,t; b(t)) \, b'(t) +
\int\limits_0^{b(t)} \frac{\partial}{\partial t} u(x, t; s) \, ds \quad (*)
$$
gives
\begin{align}
\frac{\partial}{\partial t} u(x,t) &=
\frac{\partial}{\partial t} \int\limits_0^t u(x,t;s) \, ds \\
&= \underbrace{\left. u(x,t;s) \right|_{s=t}}_0 \frac{dt}{dt} +
\int\limits_0^t \frac{\partial}{\partial t} u(x,t;s) \, ds \\
&= \int\limits_0^t u_t(x,t;s) \, ds \\
\end{align}
Applying Leibniz's rule again:
\begin{align}
\frac{\partial^2}{\partial t^2} u(x,t) &=\frac{\partial}{\partial t} \int\limits_0^t u_t(x,t;s) \, ds \\
&= \underbrace{\left. u_t(x,t;s) \right|_{s=t}}_{f(x,t)} \frac{dt}{dt} +
\int\limits_0^t \frac{\partial}{\partial t} u(x,t;s) \, ds \\
&=f(x,t) + \int\limits_0^t u_{tt}(x,t;s) \, ds
\end{align}
Derivation of equation $(*)$:
Define
$$
\varphi(x,t) := \int\limits_0^{b(t)} u(x,t;s)\; ds
$$
then using some integral mean value theorem one gets
\begin{align}
\Delta \varphi(x, t)
&= \varphi(x, t + \Delta t) - \varphi(x, t) \\
&= \int\limits_0^{b+\Delta b} u(x,t + \Delta t;s) \, ds -
\int\limits_0^b u(x,t;s) \, ds \\
&= \int\limits_0^{b} u(x,t + \Delta t;s) \, ds +
\int\limits_b^{b+\Delta b} u(x,t + \Delta t;s) \, ds -
\int\limits_0^b u(x,t;s) \, ds \\
&=\int\limits_b^{b+\Delta b} u(x,t + \Delta t;s) \, ds +
\int\limits_0^{b} \left[ u(x,t + \Delta t;s) - u(x,t;s) \right] \, ds \\
&=u(x, t + \Delta t;\sigma) \, \Delta b +
\int\limits_0^{b} \frac{u(x,t + \Delta t;s) - u(x,t;s)}{\Delta t} \, ds \, \Delta t
\end{align}
with $\sigma \in [b, b + \Delta b]$ and therefore
$$
\frac{\Delta \varphi(x, t)}{\Delta t} =
u(x, t + \Delta t;\sigma) \, \frac{\Delta b}{\Delta t} +
\int\limits_0^{b} \frac{u(x,t + \Delta t;s) - u(x,t;s)}{\Delta t} \, ds
$$
For $\Delta t \to 0$ this shrinks $\sigma \to b(t)$ and gives
$$
\varphi_t(x,t) =
u(x, t;b(t)) \, b'(t) + \int\limits_0^{b} u_t(x,t;s) \, ds
$$
Assuming $u(t,x)=T(t)X(x)$, the separated equations are
$$
\frac{T''}{c^2 T} = \lambda, \;\; \lambda = \frac{X''}{X}, \; X'(0)=X'(\pi)=0.
$$
The $X$ solutions dictate the values of $\lambda$ to be $-n^2$ for $n=1,2,3,\cdots$, and the corresponding eigenfunctions are unique up to multiplicative constants, and are given by
$$
X_n(x) = \cos(n x),\;\;\; n=0,1,2,3,\cdots.
$$
The general solution is
$$
u(x,t) = (A_0+B_0t)+\sum_{n=1}^{\infty}\left(A_n\cos(nc t)+B_n\sin(nc t)\right)\cos(n x).
$$
The constants $A_n,B_n$ are determined by the initial conditions:
$$
\cos(x) = u(x,0) = A_0+\sum_{n=1}^{\infty}A_n\cos(n x), \\
\cos^2(x) = u_{t}(x,0) = B_0+\sum_{n=1}^{\infty}nc B_n\cos(n x).
$$
The mutual orthogonality of the functions $\{ \cos(n\pi x) \}_{n=0}^{\infty}$ in $L^2[0,\pi]$ is used to determine the coefficients $A_n$ and $B_n$ in the usual manner of Fourier, which is simplified after applying the identity
$$
\cos^2(x) = \frac{1+\cos(2x)}{2}.
$$
Best Answer
First of all, you should distinguish the spacetime regions $x\le ct$ (boundary layer) and $x>ct$. When $x> ct$, the influence of the boundary is not felt, and the solution is zero due to the initial condition being zero. From now on assume $x\le ct$.
The formula $$u(x,t) = \int_{0}^{x-ct} h(-s/c)ds + G(ct-x) + G(x+ct)$$ is correct, although the integral is awkwardly written: the upper limit is negative. I would introduce $H(t)=\int_0^t h(s)\,ds$ and write $$ \int_{0}^{x-ct} h(-s/c)ds = -c \int_0^{t-x/c} h(\xi)\,d\xi = -c H(t-x/c) $$
So far, $$u(x,t) = -c H(t-x/c) + G(ct-x) + G(x+ct)\tag{1}$$
Yes. Note that plugging $t=0$ here is problematic because we work in the regime $x\le ct$. It's better to plug $x=ct$, where (by continuity) $u=0$. So, $$ 0 = G(0) + G(2x)$$ hence $G$ is identically zero.
The final answer: $$ u(x,t) = \begin{cases} -cH(t-x/c),\quad &0\le x\le ct; \\ 0,\quad & x>ct\end{cases} $$
Verification: in the boundary layer $x< ct$ we have $u_x(x,t) = h(t-x/c)$, which at $x=0$ gives $h(t)$ as required. The initial conditions are also satisfied, and the function is continuous.