I am experiencing problems when using NDSolve to find the time-dependent solution to a system of 37 coupled ODEs. ODEs contain terms which are forced through oscillatory sine functions.

This takes place in a numerical model where variables for which I want to find a time-dependent solution are the concentration of chemical species in an ocean chemistry model, while model forcing is the oscillatory evolution of the hydrologcal cycle (evaporation and river runoff water fluxes).

The code is very long and I will not copy it here, hoping that this text will be enough for some of you to guide me in a promising direction.

(1) My base setting of the NDSolve function was:

res = NDSolve[Join[ODEsys, IC], var, {t, 0, tmax}, AccuracyGoal -> 14,
   PrecisionGoal -> 14]

=> this returned the following error:

NDSolve::ndcf: Repeated convergence test failure at t == 2.4122636505104266`*^7; unable to continue.

(2) I was able to obtain a solution for simulated times tmax smaller than 10 years using the following setting:

res = NDSolve[Join[ODEsys, IC], var, {t, 0, tmax}, AccuracyGoal -> 14,
   PrecisionGoal -> 14, 
  Method -> {"EquationSimplification" -> "Residual"}, 
  MaxStepSize -> tmax/150000, MaxSteps -> Infinity]

(3) However, for longer simulation times (tmax > 10 years), I get the following error:

NDSolve::ndsz: At t == 2.9634156107275105`*^9, step size is effectively zero; singularity or stiff system suspected.

(4) Switching to the following setting which in other occasions solved stiff systems:

res = NDSolve[Join[ODEsys, IC], var, {t, 0, tmax}, AccuracyGoal -> 16,
   PrecisionGoal -> 16, WorkingPrecision -> 18, 
  Method -> {"StiffnessSwitching"}, MaxStepSize -> tmax/150000, 
  MaxSteps -> Infinity]

produces this error message:

NDSolve::precw: The precision of the differential equation ({{(TemporaryVariable$1623^\[Prime])[t]==0,<<49>>,<<24>>},{},{},{},{}}) is less than WorkingPrecision (18.`). Tensors....have incompatible shapes.

I am stuck since my knowledge of the numerous settings of the NDSolve function is very limited.

Can anyone give me a tip or direct me towards documentation where i can learn more about how to solve these kinds of issues with NDSolve ?

Thank you for having read to this point!

Vanni