If the mass ranges over 25 orders of magnitude, it is likely that other scales, such as the relevant time scale, will also vary over many orders. You can define a characteristic length scale, $y_{10}$ and frequency scale, $\omega = 2{y_{10}}\sqrt {\frac{{{C_0}}}{m}} $ and re-write your equations in dimensionless form. You will have to consider how your initial conditions scale as well.
$$m\frac{{{\partial ^2}y1\left( t \right)}}{{\partial {t^2}}} = 4{C_0}{\left( {y2\left( t \right) - y1\left( t \right)} \right)^3} - 4{C_0}y1{\left( t \right)^3}$$
$$\frac{{{\partial ^2}y2\left( t \right)}}{{\partial {t^2}}} = 4{C_0}{\left( {y1\left( t \right) - y2\left( t \right)} \right)^3} - 4{C_0}y2{\left( t \right)^3}$$
In dimensionless form, they become:
$$\frac{{{\partial ^2}y1\left( t \right)}}{{\partial {t^2}}} = {\left( {y2\left( t \right) - y1\left( t \right)} \right)^3} - y1{\left( t \right)^3}$$
$$\frac{{{\partial ^2}y2\left( t \right)}}{{\partial {t^2}}} = {\left( {y1\left( t \right) - y2\left( t \right)} \right)^3} - y2{\left( t \right)^3}$$
To estimate a characteristic length scale, $y_{10}$, I used the maximum absolute length of $y_1$ from data1 or
ymax = data1[[All, 1]] // Abs // Max (* 0.538851 *)
For the case provided, we can scale the initial condition based on the estimate of $y_{10}$ and $\omega$ as shown
$$y1'(0)=1\rightarrow y1'(0)=\frac{1}{\omega y_{10}}(Dimensionless)\Rightarrow 0.770102$$
Now, you can think of the dimensionless timescale as the number of revolutions in radians. Here is an example Mathematica code. You will need to rescale the data back to dimensional form after the NDSolve. Also note that I combined $m$ and $C_0$ into characteristic frequency $\omega$.
data2 = Block[{omega = (2 Sqrt[10] 0.5388506382979702`)/Sqrt[2] ,
y10 = 0.5388506382979702`, v10 = 1},
Reap[NDSolve[{y1''[t] == -y1[t]^3 + (y2[t] - y1[t])^3,
y2''[t] == -y2[t]^3 + (y1[t] - y2[t])^3,
y1'[0] == v10/(omega*y10), y1[0] == 0, y2'[0] == 0, y2[0] == 0,
WhenEvent[Mod[t, 2 \[Pi]] == 0,
Sow[{y1[t], y1'[t]}]]}, {}, {t, 0, 5000 (2 Pi)},
MaxSteps -> Infinity]]][[-1, 1]];
ListPlot[data3, ImageSize -> Large, PlotRange -> All,
PlotStyle -> PointSize[0.0025]]
For molecular vibrations, $\omega$ will be on the order of $10^{13-14}$ Hz, distances about $0.25x10^{-9}m$, and speed of about $25,000 m/s$ (crude estimates). With the essentially the same code and just , you can get a result quickly, although I do not claim to know the physical significance. It is more to show that you make NDSolve operate over many orders of magnitude of physical quantities by converting the problem into non-dimensional form.
data3 = Block[{omega = 10^14, y10 = 0.25 10^(-9), v10 = 25000},
Reap[NDSolve[{y1''[t] == -y1[t]^3 + (y2[t] - y1[t])^3,
y2''[t] == -y2[t]^3 + (y1[t] - y2[t])^3,
y1'[0] == v10/(omega*y10), y1[0] == 0, y2'[0] == 0, y2[0] == 0,
WhenEvent[Mod[t, 2 Pi] == 0, Sow[{y1[t], y1'[t]}]]}, {}, {t,
0, 10000 Pi}, MaxSteps -> Infinity]]][[-1, 1]];
ListPlot[data3, ImageSize -> Large, PlotRange -> All,
PlotStyle -> PointSize[0.0025]]
You can note that in data2, that the scaled ymax is indeed one. In the data3 case, they are of order 1 indicating the estimates might be in the right ball park.
