Then it must be because your theta_n has become complex. How that came to be is impossible to guess based on this information. If you use the debugger you can step up into your ODE-function and inspect its state when the error happens. A blind guess is that you have some terms (with something like square roots or logs) that should be positive but has become negative. Perhaps you can avoid this by setting the ode-options for nonnegative solutions, perhaps you have to make bigger changes to the ODE.
HTH
