There are explicit and implicit numerical solvers for ODE initial value problem.
When one solves initial value problem (or Cauchy problem) for ordinary differential equations (ODE) numerically that means we get a sequence of approximated function values y(n) starting with fist given value y(0) for a given interval and with a given discretization rate (timestep).
dy/dx = f(y, t);
y(0) = y0;
h = discretization rate: (x(n+1) - x(n))
x(end) = end of the interval
Explicit method implies that for determining next point y(n+1) we already have everything we need so we have an explicit formula for y(n+1). So for explicit method we only need to apply some formula to get next approximated point.
Implicit method works in a different way. We get an equation which references y(n+1) several times.
F(y(n+1), y(n), t) = 0
Note, that for each iteration everything except y(n+1) is already known. So the idea is to solve this equation (with independent variable y(n+1)) to find its roots. Once we have solved it, we get numerical value for y(n+1).
Explicit and implicit methods can be also one step or multi-step. One step means the formula depends only on previously calculated value y(n). Multi step method depends also on older points y(n-1), y(n-2), etc.
There is a lot of mathematical books explaining different methods, their orders (errors). I will only mention two most ancient methods here which are easy to remember and are very intuitive. However, one must know that they perform bad comparing with other advanced methods.
Explicit Euler method a.k.a. Forward Euler method
The idea is simple. We get next approximated value by using tangent at the current point.
This is an explicit method because we have explicit formula for y(n+1).
y(n+1) = y(n) + h * (f(y(n), t+1));
You can derive this method from the formula of the derivative: dy/dx = f(y(n), t+1) =
(y(n+1) – y(n)) / h
Geometrically it works like this:

Implicit Euler method a.k.a Backward Euler method
y(n+1) = y(n) + h * (f(y(n+1), t+1));
Note on both side of equation we have y(n+1) which is not known.
So we need to make reform this equation so it would look like
F(y(n+1), y(n), t) = 0
We get
y(n+1) - y(n) - h * (f(y(n+1), t+1)) = 0;
Here we are really depend on the given function f because sometimes it is possible to make an explicit formula for y(n+1) but most of the times it is not possible.
So, imagine you are finding roots of an equation (substitute y(n+1) with x):
x - C1 - C2 * (f(x, C3)) = 0;
where x is independent variable and all other constants are known. Depending on the complexity of f we can get some analytical solution, but for the general case this equation is solved numerically. For the roots finding algorithm one can use simple bisection method but Newton method is more appreciated since it gives better results much faster.
Now we get approximated root of the equation so we know y(n+1) which is our next approximated point.
Geometrically backward Euler is also very intuitive:

Both Euler methods mentioned here have O(h) error order. That means if we want to make error twice better we need to make timestep twice smaller. For example a Runge-Kutta method of order O(h^4) allows to have 16 times less error when we divide timestep in half.
Now about the implementation of these methods.
There are abbreviations called PECE or PECECE or maybe you can find a PECECECE.
PE means there is a predicted value (expectation).
CE means there is a corrected value (estimation).
PECE also means predictor corrector scheme.
This is a way of using implicit methods like backward Euler or Adams-Moulton.
Instead of solving equation for y(n+1), we use explicit method of the same order to calculate predicted value.
Then, we insert y(n+1) into the implicit formula and we get y(n+1) which is corrected.
This idea reduces computational costs since otherwise for implicit methods you need to solve arbitrary equation (linear of non linear) with numerical algorithm.
I have implemented in Matlab explicit and implicit Euler, Adams-Moulton, Heun and Runge-Kutta of order 4 which you can download all in one archive as ode_matlab.zip. There is also Newton root finding algorithm implemented in this archive since for implicit method we need to solve equation numerically to get next approximated value.
Explicit methods are computationally faster since they require less function evaluations. Implicit methods are slower but more stable again stiff equations. That means on some problems explicit method may require very small discretization rate (timestep) where implicit method will produce appropriate solution with much greater timesteps.
Thanks for the code. Do you have a C++ version?