The shape evolution and shape coexistence phenomena in neutron-rich nuclei at $Napprox60$, including Kr, Sr, Zr, and Mo isotopes, are studied in the covariant density functional theory (DFT) with the new parameter set PC-PK1. Pairing correlations are treated using the BCS approximation with a separable pairing force. Sharp rising in the charge radii of Sr and Zr isotopes at N=60 is observed and shown to be related to the rapid changing in nuclear shapes. The shape evolution is moderate in neighboring Kr and Mo isotopes. Similar as the results of previous Hartree-Fock-Bogogliubov (HFB) calculations with the Gogny force, triaxiality is observed in Mo isotopes and shown to be essential to reproduce quantitatively the corresponding charge radii. In addition, the coexistence of prolate and oblate shapes is found in both $^{98}$Sr and $^{100}$Zr. The observed oblate and prolate minima are related to the low single-particle energy level density around the Fermi surfaces of neutron and proton respectively. Furthermore, the 5-dimensional (5D) collective Hamiltonian determined by the calculations of the PC-PK1 energy functional is solved for $^{98}$Sr and $^{100}$Zr. The resultant excitation energy of $0^+_2$ state and E0 transition strength $rho^2(E0;0^+_2rightarrow0^+_1)$ are in rather good agreement with the data. It is found that the lower barrier height separating the two competing minima along the $gamma$ deformation in $^{100}$Zr gives rise to the larger $rho^2(E0;0^+_2rightarrow0^+_1)$ than that in $^{98}$Sr.