To update you all... I am even more baffled, because the solver seems to correctly find solutions to the system if I set
ft=[0:0.01:5e4];
f=0.01*cos(0.002*ft);
This is a massive simplification of the forcing vector I would like to use, but I used it as a test to see if my ode code was correct... looks like it is, but the problem may lie with the vector I am trying to force the system with.
Perhaps my vector has a timestep that is too small (0.025ms), but it works for a cosine vector of the same timestep.
So I am thinking that it is the forcing function that is the problem and leading to the error:
Warning: Matrix is singular, close to singular or badly scaled.
Results may be inaccurate. RCOND = NaN.
But the system, when simulated in XPPAUT, is able to find a solution...