Second answer -- OK, the first answer was hogwash (the curious can inspect the edit history). It pays sometimes to write out the equation and think about it first. The code for this one looks so complicated, but the equations basically have the form (here d = 80000)
rc'[t] == 2.05594*10^-10/rc[t] - 8189.14 rc[t] + 80000. rm[t],
rm'[t] == -80000. rc[t] + 2.05594*10^-10/rm[t] + 2189.14 rm[t]
We see that there are singularities along both axes rc = 0 and rm = 0. So the solution should not be crossing the axes at all. Here is a simplified vector field that will illustrate the problem NDSolve faces.
StreamPlot[{1/rc - rc + 2 rm, -2 rc + 1/rm + rm},
{rc, -3, 3}, {rm, -4.5, 4.5}, Axes -> True]
As a trajectory moves clockwise away from an axis, it approaches the next axis before turning toward an equilibrium point or out toward infinity, depending on the quadrant. The turn can be rather sudden. If NDSolve, taking a discrete step, happens to pass over the axis far enough, the solution will continue in the same direction. This will lead to a computed solution that is oscillatory. With the OP's coefficients, the regions near the axes where the phase vectors turn are very narrow and the OP's NDSolve code often jumps over it. When it doesn't jump over an axis, the oscillations are interrupted, until it does again. In short it's a step-size issue.
What should happen to a solution that starts in the first quadrant is that it should converge on the equilibrium position.
Solution
Here is the set-up of the problem. I put a variable prec to specify the precision, since it took some experimentation to find a robust setting. It may be that prec needs to be adjusted further to handle a particular use-case.
(*Solve for dynamics*)
d = 80000;
prec = 40;
de := SetPrecision[
{rc'[t] == rate*B/w1[V]^2/2/rc[t] - rate*wc[V]*A*rc[t]/w1[V]/2 + d*rm[t],
rm'[t] == rate*B/w1[V]^2/2/rm[t] + rate*wm[V]*A*rm[t]/w1[V]/2 - d*rc[t],
rm[0] == rm0, rc[0] == rc0} /. {B -> 3.4*^-4, A -> 3*^-3},
prec]
The equilibria:
NSolve[DeleteCases[de, _?(FreeQ[#, t] &)] /. {rc'[t] -> 0., rm'[t] -> 0.},
{rc[t], rm[t]}, Reals]
(*
{{rc[t] -> -2.24157*10^-7, rm[t] -> -1.14809*10^-8},
{rc[t] -> 2.24157*10^-7, rm[t] -> 1.14809*10^-8}}
*)
There are several strategies and parameters of NDSolve that figure into obtaining an accurate solution. Some of them required experimentation, and I left some of this experimental code in, just in case someone wants to play with it. It is difficult to describe the interdependence of the elements of the code, so I will explain each one separately.
WhenEvent -- This can be used to detect either rc or rm changing sign at a step taken by NDSolve. If this happens, NDSolve should back up and start over. The action "RestartIntegration" causes this.
"LocationMethod" -- The setting "LocationMethod" -> "StepBegin" causes NDSolve to restart the integration from the beginning of the current step. If we restarted at the zero crossing, we would get a divide by zero error if using arbitrary precision. If we were using MachinePrecision, we would get, a more or less randomly, either the solution staying on the correct side of the axis or skipping across.
"StartingStepSize" -- When NDSolve starts, this is irrelevant. Its purpose is to reset the step size after the WhenEvent action of "RestartIntegration". (This is the clever bit, IMO, if there is any cleverness at all here.)
"IntegrateEvent" -- This makes NDSolve take smaller steps near an event. This seemed to be important sometimes.
WorkingPrecision -- This problem seems to have trouble with round-off error and numbers rather close to zero. It is important to have a high working precision, so that the differential equation and events are calculated accurately.
AccuracyGoal -- It is important that this be high. The coordinates of the equilibrium are around 10^-8 and sometimes smaller with other parameter settings. In some use-cases this may need to be raised.
PrecisionGoal -- By default, PrecisionGoal is set to be half of WorkingPrecision. When a computation is numerically challenging, it is sometimes important to increase the ratio of WorkingPrecision to PrecisionGoal. This seems to be the case here. On the other hand, the higher the PrecisionGoal is, the more work NDSolve will do (i.e., it will take longer and use more memory). While I used round numbers, AccuracyGoal exceeds PrecisionalGoal by about how many leading zeros there are in smallest coordinate of the equilibrium position.
MaxStepSize -- This particular problem seems hard to control with MaxStepSize but I left it in, commented out, in case it might be useful in a particular use-case.
Reap/Sow -- This is also part of the code for experimenting. StepMonitor can be used to see what steps are taken in computing the solution. One can inspect the values of t, rc[t] and rm[t] actually used by NDSolve. I used it to monitor the step size. Uncomment the GridLines options in the Plot to see the steps.
Solution:
{mysolve, tvals} = Reap[NDSolve[{
de,
WhenEvent[rc[t] < 0, "RestartIntegration",
"LocationMethod" -> "StepBegin", "IntegrateEvent" -> True],
WhenEvent[rm[t] < 0, "RestartIntegration",
"LocationMethod" -> "StepBegin", "IntegrateEvent" -> True]},
{rc, rm}, {t, 0., 0.04},
MaxSteps -> Infinity,
StartingStepSize -> 1*^-11, (*MaxStepSize->1*^-6,*)
AccuracyGoal -> 20, PrecisionGoal -> 10, WorkingPrecision -> prec,
StepMonitor :> Sow[t]]];
Plot:
(* to determine the plot range *)
{cPR, mPR} = {Max[Abs@rc[t] /. #], Max[Abs@rm[t] /. #]} &@
NSolve[DeleteCases[de, _?(FreeQ[#, t] &)] /. {rc'[t] -> 0., rm'[t] -> 0.},
{rc[t], rm[t]}, Reals];
(*Plot dynamics*)
Plot[Evaluate[rc[t] /. mysolve], {t, 0.0, 0.04},
PlotRange -> {-0.1 cPR, 1.5 cPR}, AspectRatio -> 0.25,
ImageSize -> Full (*,GridLines->{First@tvals,None}*)]
Plot[Evaluate[rm[t] /. mysolve], {t, 0.0, 0.04},
PlotRange -> {-0.1 mPR, 1.5 mPR}, AspectRatio -> 0.25,
ImageSize -> Full]
