53 std::vector<double>& rCurrentYValues,
54 std::vector<double>& rNextYValues)
61 const double delta = 1.0e-8;
65 if (
mEvalF.size() != num_equations)
67 mEvalF.resize(num_equations);
69 mTemp.resize(num_equations);
70 mYinit.resize(num_equations);
73 rNextYValues = rCurrentYValues;
79 for (
unsigned i=0; i<num_equations; i++)
81 rNextYValues[i]=rNextYValues[i]+delta;
84 rNextYValues[i]=rNextYValues[i]-delta;
87 for (
unsigned i=0; i<num_equations; i++)
91 rNextYValues[i]=rNextYValues[i]+0.5*
mEvalF[i]*timeStep;
99 for (
unsigned i=0; i<num_equations; i++)
101 ysave = rNextYValues[i];
102 rNextYValues[i]=
mYinit[i];
106 rNextYValues[i]=rNextYValues[i]+delta;
109 rNextYValues[i]=ysave;
113 for (
unsigned i=0; i<num_equations; i++)