The solution is OK, if you know that you have to perform a time propagation of exactly \$t\$ with that same Hamiltonian and you have to do so countless times.
As someone else has suggested, ditch all the matmul and conjg and use only BLAS calls.
If you apply the propagator to a small number of vectors, however, computing the full propagator is not efficient, due to the matrix-matrix multiplication it entails. Instead, you should evaluate its action on the wavefunction you are propagating by means of three matrix-vector multiplications:
$$e^{-i\mathbf{H}dt}\mathbf{c} = \mathbf{U} e^{-i \mathbf{E}dt}\mathbf{U}^\dagger \mathbf{c}$$
that is
$$\mathbf{c}_1 = \mathbf{U}^\dagger \mathbf{c},\qquad \mathbf{c}_2 = e^{-i \mathbf{E}dt} \mathbf{c}_1,\qquad \mathbf{c} = \mathbf{U} \mathbf{c}_2$$
If you need to propagate multiple different times, of course, diagonalizing the Hamiltonian for every time is definitely not a good idea. Instead, you should do the diagonalization only once, and follow the triple-multiplication procedure above, whose computational cost scales as \$N^2\$ rather than \$N^3\$.
If your hamiltonian were to ever become time-dependent, \$\mathbf{H}=\mathbf{H}(t)\$, you must break the propagation in steps
$$\mathbf{c}(t+dt) = \mathbf{U}(t+dt,t)\mathbf{c}(t),$$
where \$\hat{U}(t+dt,t)\$ is the time-ordered exponential propagator,
$$\mathbf{U}(t+dt,t)=\hat{T}e^{-i\int d\tau\mathbf{H}(\tau)}$$.
An excellent unitary second-order propagator is the mid-point exponential
$$ c(t+dt) = e^{-i \mathbf{H}(t+dt/2) dt} c(t).$$
To evaluate this step, there are many ways that are much better than a brutal direct diagonalization at each time step. First, if you can afford the diagonalization of the full Hamiltonian and if the Hamiltonian has the form
$$\mathbf{H}(t) = \mathbf{H}_0 + F(t) \mathbf{H}_I$$
where both \$\mathbf{H}_0\$ and \$\mathbf{H}_I\$ are time-independent Hermitean operators and \$F(t)\$ is a real function, then the best approach is to split symmetrically the propagator:
$$e^{-i \mathbf{H}(t+dt/2) dt} = e^{-i \mathbf{H}_0 dt/2}e^{-i \mathbf{H}_I F(t+dt/2) dt}e^{-i \mathbf{H}_0 dt/2} + o(dt^3),$$
where \$ o(dt^3)\$ is the difference between the two approximations to the time-step propagator. The mid-point propagator was already accurate only to second order anyway, so you are not losing any accuracy. With this arrangement, you need to perform the diagonalization only once, at the very beginning of the program, separately for \$\mathbf{H}_0\$ and \$\mathbf{H}_I\$,
$$\mathbf{H}_0 = \mathbf{U}_0 \mathbf{E}^{(0)} \mathbf{U}_0^\dagger,\qquad\mathbf{H}_I = \mathbf{U}_I \boldsymbol{\Delta} \mathbf{U}_I^\dagger,$$
where \$\mathbf{E}^{(0)}_{ij}=E^{(0)}_i\delta_{ij}\$, \$\boldsymbol{\Delta}_{ij}=\Delta_i\delta_{ij}\$ are diagonal matrices.
In this way, your propagator becomes
$$U(t+dt,t) = \mathbf{U}_0 e^{-i \mathbf{E}^{(0)} dt/2}\mathbf{U}_0^\dagger\mathbf{U}_I e^{-i \boldsymbol{\Delta}_I F(t+dt/2) dt}\mathbf{U}_I^\dagger\mathbf{U}_0 e^{-i \mathbf{H}_0 dt/2}\mathbf{U}_0^\dagger + o(dt^3),$$
You can pre-compute the matrix \$\mathbf{O}=\mathbf{U}_0^\dagger\mathbf{U}_I\$, of course, in which case
$$U(t+dt,t) = \mathbf{U}_0 e^{-i \mathbf{E}_0 dt/2}\mathbf{O} e^{-i \mathbf{E}_I F(t+dt/2) dt}\mathbf{O}^\dagger e^{-i \mathbf{H}_0 dt/2}\mathbf{U}_0^\dagger + o(dt^3).$$
If your time step is constant, you can neglect the first half free propagation, and merge the second with the first of the following step, thus obtaining
$$U(t+dt,t) \simeq \mathbf{U}_0 e^{-i \mathbf{E}_0 dt}\mathbf{O} e^{-i \mathbf{E}_I F(t+dt/2) dt}\mathbf{U}_I^\dagger.$$
In this way, you have replaced a diagonalization (order \$N^3\$) with three consecutive matrix-vector multiplications (order \$N^2\$).
If applicable to your case, therefore, this change will result in a wild speed up.
If either \$H_0\$ or \$H_I\$ are too large to be diagonalized (they are not equivalent, since \$H_0\$ is typically block diagonal and hence it is sufficient to diagonalize its diagonal blocks), then you can still use the symmetric splitting as a way to precondition the matrix, and use an iterative method (look up Krylov spaces) to determine the action of the operator on your target function, rather than the opearator over the whole Hilbert space, which is almost always an overkill.
PS: For folks not familiar with the simple syntax of modern Fortran, I suggest quick modern Fortran tutorial.
