Классичкский метод Рунге-Кутта 4-го порядка (RK4)

Добрый день)

Понадобился мне тут на днях класс, для решения систем обыкновенных диф. уравнений (ОДЕ) методом, желательно наиболее устойчивым)) Погуглив на эту тему, было решено остановиться на классическом методе Рунге-Кутта 4-го порядка.
В общем выкладываю реализацию, которая первая пришла в голову.. Но мне что то она не оч нравиться((
Как возможно сделать наиболее правильно, чтоб с одной стороны этим классом было удобно и интуитивно понятно пользоваться, а с другой чтоб производительность не страдала, а напротив?)

abstractode.h

Код

C++ (Qt)
#ifndef ABSTRACTODE_H
#define ABSTRACTODE_H
 
#include <vector>
 
class AbstractODE
{
public:
    AbstractODE();
    virtual ~AbstractODE();
    std::vector<double> operator()(const std::vector<double> &initPoint,
                                   const double &t0, const double &t, const double &h=0.001) const;
    virtual int dimension() const = 0;
 
protected:
    virtual std::vector<double> function(const std::vector<double> &, const double &) const = 0;
 
private:
};
 
#endif // ABSTRACTODE_H

abstractode.cpp

Код

C++ (Qt)
#include "abstractode.h"
#include <cmath>
 
AbstractODE::AbstractODE()
{
}
//-----------------------------------------------------------------------------
 
AbstractODE::~AbstractODE()
{
}
//-----------------------------------------------------------------------------
 
std::vector<double> AbstractODE::operator ()(const std::vector<double> &initPoint,
                                             const double &t0, const double &t, const double &h) const
{
    std::vector<double> x = initPoint;
    if ((int)x.size() != dimension())
        x.resize(dimension());
 
    std::vector<double> k1(dimension());
    std::vector<double> k2(dimension());
    std::vector<double> k3(dimension());
    std::vector<double> k4(dimension());
 
    std::vector<double> cx(dimension());
 
    double _t = t0;
    double _h = fabs(h);
    if (t < t0)
        _h *= -1;
 
    bool isStop = false;
 
    do {
 
        if (fabs(t - _t) <= fabs(h)) {
            _h = t - _t;
            isStop = true;
        }
 
        for (int i = 0; i < dimension(); i++) {
            k1[i] = function(x, _t)[i];
        }
        for (int i = 0; i < dimension(); i++) {
            cx[i] = x[i] + 0.5 * _h * k1[i];
        }
        double tau = _t + 0.5 * _h;
        for (int i = 0; i < dimension(); i++) {
            k2[i] = function(cx, tau)[i];
        }
        for (int i = 0; i < dimension(); i++) {
            cx[i] = x[i] + 0.5 * _h * k2[i];
        }
        for (int i = 0; i < dimension(); i++) {
            k3[i] = function(cx, tau)[i];
        }
        for (int i = 0; i < dimension(); i++) {
            cx[i] = x[i] + _h * k3[i];
        }
        tau = _t + _h;
        for (int i = 0; i < dimension(); i++) {
            k4[i] = function(cx, tau)[i];
        }
 
        for (int i = 0; i < dimension(); i++) {
            x[i] = x[i] + _h * (k1[i] + 2.0 * (k2[i] + k3[i]) + k4[i])/6.0;
        }
 
        _t += _h;
 
    } while (!isStop);
 
    return x;
}
//-----------------------------------------------------------------------------

Используется енто так:
Вначале наследуемся от AbstractODE и переопределяем две виртуальные функции:
1) int dimension() const;
2) std::vector<double> function(const std::vector<double> &x, const double &) const;
пример:

Код

C++ (Qt)
#include <iostream>
#include <fstream>
#include "abstractode.h"
#include <cmath>
 
using namespace std;
 
class MyODE : public AbstractODE
{
public:
    int dimension() const { return 2; }
protected:
    std::vector<double> function(const std::vector<double> &x, const double &) const;
 
};
 
inline std::vector<double> MyODE::function(const std::vector<double> &x,
                                           const double &/*t*/) const
{
    static std::vector<double> res(dimension());
 
    res[0] = x[1];
    res[1] = -x[0];
 
    return res;
}
 
 
int main(int argc, char *argv[])
{
    MyODE myODE;
    std::vector<double> initPoint(myODE.dimension());
    initPoint[0] = 0.5;
    initPoint[1] = 0.0;
 
    const double t0 = 0.0;
    const double t = 10.0;
 
    std::vector<double> x = myODE(initPoint, t0, t, 1e-2);
 
    cout << " x[0] = " << x[0] << endl;
    cout << " x[1] = " << x[1] << endl;
 
    return 0;
}

Вопрос: Как улучшить?

Заранее благодарю)