I believe that the best way to solve this is to draw the $x$-$t$ plane and use the characteristic lines and the fact that the solution to this equation is constant along the characteristic lines.
The first section of the $x$-$t$ plane is when $x>t$ (right of the line with slope 1), which corresponds to information from the initial condition, so $u=0$ in this region.
The second region is where information is carried from the BC $u=1$ from $t\in(0,1)$. In this region, the characteristics have slope $1/3$ in the $x$-$t$ plane, so the characteristics will collide with $x=t$, forming a shock. Nevertheless, we have $u=1$ when $x<t$ and $x>3(t-1)$.
For the third region, which corresponds to $x<3(t-1)$, the slopes of the characteristics vary depending on where they start on the $t$-axis, the value of which we will call $s$. The characteristic lines then have the form $t=x/(1+2s)+s$ and along this characteristic, $u=s$. We can solve this first equation for $s$, which yields a quadratic equation with 2 real roots. One of them is not feasible (remember we need $s>1$) and the other is, so we can write $u(x,t)=s=(-1 + 2 t + \sqrt{(2t+1)^2 - 8 x})/4$
$$
u(x,t) =
\begin{cases}
0, & x\geq t \\
1, & x<t \text{ and }x\geq 3(t-1) \\
\frac{-1 + 2 t + \sqrt{(2t+1)^2 - 8 x}}{4}, & x<3(t-1)
\end{cases}, \ \forall x,t>0.
$$
It seems I was confused or made mistakes the first time trying to solve this, so I will post my solution now that it makes sense.
Assume $u(x,t)=F(x+2t)+G(x-2t)$. Then, the initial conditions give us:
$$
u(x,0)=F(x)+G(x)=e^{-x}\\
u_t(x,0)=2F'(x)-2G'(x)=2e^{-x}
$$
which hold for all $x>0$.
Hence, in the last equation, we may divide both sides by 2 and integrate with respect to x, then solve the resulting system:
$$
\begin{cases}
F(x)+G(x)=e^{-x}\\
F(x)-G(x)=-e^{-x}+C
\end{cases}
$$
where $C$ is some constant.
We thus obtain $F(x)=C/2$ and $G(x)=e^{-x}-C/2$. Since the only condition is that $x$ is positive, we may replace it by any positive quantity (I think this is where I was previously confused). Thus, we obtain $u(x,t)=F(x+2t)+G(x-2t)=C/2+e^{-(x-2t)}-C/2=e^{2t-x}$ which holds for all $x>2t$. It is clear, as I noted in my original post, that this corresponds to the same solution that we get from d'Alembert's formula.
Now we handle the case when the argument to $G$ is negative. For this, we need to use the boundary condition.
We have $u_x(0,t)=F'(2t)+G'(-2t)=-\cos t$ for all $t>0$. Make the substitution $z=-2t$ to obtain $G'(z)=-\cos(z/2)-F'(-z)=-\cos(z/2)$ by noting that $F'(-z)=0$ from our previous work. Integrating, we obtain $G(z)=-2\sin(z/2)+\tilde{C}$.
Applying the continuity condition, we must have $G(0)=\tilde{C}=1-C/2$. Thus, we have $u(x,t)=F(x+2t)+G(x-2t)=1-2\sin(x/2-t)$ for $0<x<2t$.
Hence, the complete solution is:
$$u(x,t)=\begin{cases}1-2\sin(x/2-t),&0<x\leq 2t\\e^{2t-x},&x>2t\end{cases}.$$
We can now verify that the initial conditions, boundary condition, and continuity are satisfied.
Best Answer
I'm finding your last comment a little hard to understand, so I will just demonstrate the solution below. We have the initial data $$u(x,0)=u_0(x)=(|x|-\pi)^2 \\ (\partial_tu)(x,0)=\dot{u}_0(x)=0$$ So, $$u(x,t)=\frac{1}{2}\big(u_0(x-ct)+u_0(x+ct)\big)+\frac{1}{2c}\int\limits_{x-ct}^{x+ct}\dot{u}_0(s)\mathrm ds \\ u(x,t)=\frac{1}{2}\big((|x-ct|-\pi)^2+(|x+ct|-\pi)^2\big)$$ Because integrating zero gives you zero. One can find that this solution actually unbounded, in other words it is non-physical. In order to get a physical solution, one would need to drop the initial condition $(\partial_tu)(x,0)=0$ and replace it with the boundary conditions $$u(-\pi,t)=u(\pi,t)=0$$ And this would require a different solution method, most likely using S.O.V and Sturm-Liouville theory.