Let $A : [0,T]\to\mathbb R^{n\times n}$ be a matrix function such that each $A(t)$ is normal, i.e., $A(t)^TA(t) = A(t)A(t)^T$ for all $t\in [0,T]$. Now, consider the ODE system $\dot x(t) = A(t)x(t)$. We know that there exists a $n\times n$ fundamental matrix function $\Phi(t)$ such that the columns of $\Phi$ are linearly independent solutions of the ODE, that is, $\dot\Phi(t) = A(t)\Phi(t)$ for $t\in [0,T]$ and fixing the value $\Phi(0)$ also fixes the fundamental system.
My question:
If $A$ is normal as above, does the same hold for the solution of $\dot\Phi(t) = A(t)\Phi(t)$ with $\Phi(0) = I_n$?
For constant $A$ this is true as the solution is given by $\Phi(t) = e^{tA}$, which is easily seen to be normal for each $t$. It's also true when $A(t)$ is diagonal for each $t$.
No, this is not true. Let $$A(t)=\begin{cases}\begin{bmatrix}0&\log{(2)}\\\log{(2)}&0\end{bmatrix} & 0\leq t<1 \\ \begin{bmatrix}0&\frac{\pi}{4}\\-\frac{\pi}{4}&0\end{bmatrix} & t\geq1\end{cases}$$ Then $\Phi(0)=I_2$, $\Phi(1)=\frac{1}{4}\begin{bmatrix}5&3\\3&5\end{bmatrix}$ (so far so good) but $\Phi(2)=2\sqrt{2}\begin{bmatrix}4&4\\-1&1\end{bmatrix}$, which is not normal.
To see how I got this, note that this sort of phenomenon is generic: assuming $A$ constant, on a finite interval, the solution $\Phi$ is continuous in $A$. So if we can approximate $A$ in the uniform norm, then $\Phi$ will depend continuous on the approximation (in, e.g., the uniform norm on $(\Phi,\dot{\Phi})$). In particular, we can approximate $A$ by a piecewise-constant function.
So suppose $A$ attains value $M_1$ on $[0,t_1)$ and $M_2$ on $[t_1,t_2)$, etc. Then $$\Phi(t_2)=e^{M_2t_2}e^{M_1t_1}\Phi(0)$$ Now, if $A$ is very nice, and $M_1$ and $M_2$ commute, then we can reduce this to $e^{M_2t_2+M_1t_1}\Phi(0)$. If they don't commute, this value is off...but not by much; c.f. the Baker-Campbell-Hausdorff formula. In any case, this is a real danger: a sum of normal matrices need not be normal! A classic example is $$\begin{bmatrix}0&a+b\\a-b&0\end{bmatrix}=\begin{bmatrix}0&a\\a&0\end{bmatrix}+\begin{bmatrix}0&b\\-b&0\end{bmatrix}$$ The remainder is just messing around with the constants.
Edit: the OP has asked for a $C^1$ example. The key idea is that the solution varies continuously in $A$, as measured in the $L^1$ norm of $A$. More precisely: let $B$ match $A$ on $[0,1-\epsilon]$ and $[1+\epsilon,2]$, and quadratically interpolate between. Say $\Phi_A$ and $\Phi_B$ solve the obvious systems; clearly $\|A-B\|$, $\|\Phi_A\|$, and $\|\Phi_B\|$ are all bounded (in norm) on $[0,1+\epsilon]$ by, say, $M$. Then $$\|\Phi_B(1+\epsilon)-\Phi_A(1+\epsilon)\|\leq M((1+\epsilon)-(1-\epsilon))\sup_{s\in(1-\epsilon,1+\epsilon)}{\|A(s)-B(s)\|}\leq2M\epsilon$$ On $[1+\epsilon,2]$, we have again $A=B$, so that $$\frac{d}{dt}(\Phi_A-\Phi_B)=A(\Phi_A-\Phi_B)$$ Taking norms and solving, on $[1+\epsilon,2]$, that $2M\epsilon$-sized difference grows by a constant factor of $e^{\|A\|(1-\epsilon)}$. But $$\lim_{\epsilon\to0^+}{2M\epsilon e^{\|A\|(1-\epsilon)}}\to0$$ and there are no normal matrices near $\left[\begin{smallmatrix}4&4\\-1&1\end{smallmatrix}\right]$. (To see that last claim, note that $XX^{\mathsf{T}}-X^{\mathsf{T}}X$ is continuous in $X$.)