C++ version of DUMKA3 is a template function dumka3.cpp.
Idea to use templates was taken from IML++ library. In order to use
templates you need to have classes VectorFunction and Vector

VectorFunction  implements VectorFunctionBase  interface with
three virtual functions, which have to be implemented

void compute(Vector& y,
double t, Vector& result)  calculates right hand side of the ODE
y'=f(y,t)

double cour(Vector& y,
double t)  used if "isComputeEigenValue=false" to

calculate
h_{euler}  step which you would use in explicit Euler method.

If
"isComputeEigenValue=true" internal eigenvalue
function will be used to evaluate h_{euler}=2/max(eigenvalue).

Element GetNonZeroElement()
 can be used if "isComputeEigenValue=true" in eigenvalue
function, if

iterations
of power method have to be restarted, but y,f(y) equals to zero,

in this
case program need to initialize random nonzero vector, because

vector's
Element is unknown in advance

for
example if Vector contains "double"(s) you can return
rand()%100/100.

but if
Vector contains "Point"(s) you can

return
Point(rand()%100/100.,rand()%100/100.,rand()%100/100.)


template <class Vector,
class Element>

class VectorFunctionBase {

public:

virtual void
compute(Vector& y, double t, Vector& result) = 0;

virtual double
cour(Vector& y, double t) = 0;

virtual Element
GetNonZeroElement() = 0;

};

Vector class should have function size() and operator[]
(for example class vector<double> has these functions).

Each Element of the Vector should have operator* (for
example if your Vector has elements "double" everything is
OK, but if you use "Point" you should define operator*)
This is the first C++ version of dumka3, and as everything new it will
be changed and modified in a future, and your
suggestions/complains/questions would be very welcome.
Both fortran version of dumka3.for and C++
version of dumka3.cpp require additional
function which determine maximal step of explicit Euler method. It can
be evaluated by value cou=2/M, where M is upper boundary of maximal
eigen value. For example if you solve ODE y'=5*y, then take
cou=2/5 or something smaller, for example 2/5.1. You can use also norm
of jacobian of the right hand side function.