The effect of the rare-gas matrices on the vibrational spectra of solvated molecular ions was investigated. An approach based on the full-dimensional coupled-clusters description of the gas phase ion augmented with the Monte Carlo averaging of frequency shift due to the interaction with the matrix is presented. The calculated nu(R)(CO) of the Cu+CO ion in the Ne matrix at 2229 cm(-1) agrees within 5 cm(-1) with the experimentally determined value. The remaining small discrepancy between theory and experiment is due to the error in harmonic frequency calculations of the gas phase Cu+CO ion at the coupled-clusters level of theory. Calculated frequency shift due to the effect of Ne matrix allows us to predict the gas-phase harmonic omega=2272 cm(-1) and anharmonic nu=2243 cm(-1) vibrational frequencies for Cu+CO. A large difference (-31 cm(-1)) in Deltanu(R)(CO) stretching frequencies between Ne and Ar matrices was found. The Deltanu(R)(CO) does not depend on the cavity size and on the temperature. Major contribution to Deltanu(R)(CO) comes from five Ar atoms. Together with the CO ligand these five atoms are octahedrally coordinated to the Cu+ ion. The effect of the basis set superposition error was found to be important for the calculations of Deltanu(R)(CO) in the rare-gas matrix.