ode — Создание библиотеки классов AbstractOdeSolver для численного решения дифференциальных уравнений (c ++)

Я пытаюсь сделать и од решатель, используя с ++ и численные методы (Эйлер, Хеун и Рунге Кутта). Сначала я создал и абстрагировал класс для требований решения задач, затем я создал отдельный класс для каждого решателя, унаследованного от ABClass. с кодом нет проблем, за исключением того, что он работает только с дифференциальными уравнениями первого порядка при dydt = y
однако мне нужно расширить его, чтобы иметь возможность решать систему дифференциальных уравнений первого порядка как dydt = x-y & dxdt = y вместе

вот файл заголовка

    #ifndef ABSTRACTODESOLVER_H_INCLUDED
#define ABSTRACTODESOLVER_H_INCLUDED
#include <iostream>

using namespace std;

class AbstractOdeSolver
{
private:
double stepsize;
double initialTime;
double finalTime;
double initialValue;
public:

double (*RHS)(double, double); //pointer to function

void SetStepSize(double h);
void SetIntervalTime(double t0, double t1);
void SetInitialValue(double y0);
void setRHS(double (*pRHS)(double, double)); //set pointer function

double getStepSize(){return stepsize;}
double getInitialTime(){return initialTime;}
double getFinalTime(){return finalTime;}
double getInitialValue(){return initialValue;}

virtual void SolveEquation() = 0;
};

class EulerSolver : public AbstractOdeSolver
{
public:
virtual void SolveEquation();
};

class HeunSolver : public AbstractOdeSolver
{
public:
virtual void SolveEquation();
};

class RungeKuttaSolver : public AbstractOdeSolver
{
public:
virtual void SolveEquation();
};
#endif // ABSTRACTODESOLVER_H_INCLUDED

и это исходный код для одного решателя:

    void EulerSolver::SolveEquation(){

double yNew = 0.0;
double yOld = getInitialValue(); // y0 initial value of Y
double tInit = getInitialTime(); // t0 initial value of T
double tFinal = getFinalTime();
double h = getStepSize();

for (double i = tInit; i <= tFinal; i += h){

yNew = yOld + (h * RHS(tInit,yOld));
yOld = yNew;
tInit += h;
cout << left << setw(5) << tInit << " " << setw(5) << yNew << endl;
}

}

и это программа:

    double Func(double t, double y){
return y;
}

EulerSolver euler1;
euler1.SetIntervalTime(0,0.6);
euler1.SetInitialValue(1);
euler1.SetStepSize(0.1);
euler1.setRHS(Func);
euler1.SolveEquation();

здесь я использую указатель функции RHS, чтобы позволить решателю использовать функцию каждый раз в итерации, однако мне нужно иметь более одной функции для решения уравнений более высокой степени, и проблема, которую я не знаю, как сделать два отдельных функции со связанными переменными (dydt = xy & dxdt = y) и затем наличие указателя на каждый из них и возвращение ответа в массиве!
Любые идеи будут отличными.