Thanks for all your comments. They were really helpful.
I. It turned out that the problem is not with infinities or nan values. They never occured during the calculation.
However, I always got this message: Failure at t=4.164915e-10. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.479675e-24) at time t
before the shortening of vector 't'. My time values were in the order of 1e-12 seconds. I thought this is small value so it would better if I reformulate my ODE system such that it uses picoseconds instead of 1e-12 seconds. After the implementation of this change I have not got any failure.
But I got very long running time. Even after 6 hours the ode solver (I tried all) did not even calculated the half of the points I needed.
As you said before, the ode solver decides the stepsize during the calculation, and there I found a problem: line 542-543 in ode15s :
% Predict a solution at t+h.
tnew = t + h;
The h value went below 1e-12 while the solver calculated a point between -3 and -4 ps. But I do not need better resolution than 1 ps unit.
It seems that some sets of parameters (I do not mean here the initial paramteres of the ODE system) gives the solver very hard time to make reasonable stepsize. This part I still do not understand.
II. The solution: using a solver with the fixed stepsize. I chose ode5 and it worked remarkably well. Not only for solving the ODE system but also using within the jacobianest function.
I performed the fitting and I got a consistent story in the end. Thanks again very much, I have really learnt a lot.