Gear Method.
Let's examine the following SLDE:

y₁′ = y₂
y₂′ = y₁
y₁(0)=0  and  y₂(0)=1

NB: we expect sin(x)/cos(x) as the result.

#include <stdio.h>
#include <math.h>
#include "libnum/gear.h"

using namespace NumericMethods;

bool derfun1(int n, double x, const double y[], double dy[])
{
  dy[0] = y[1];
  dy[1] = -y[0];
  return true;
}

bool jacobian1(int n, double x, const double y[], double j[])
{
  j[0] = 0; j[1] = 1; j[2] = -1; j[3] = 0;
  return true;
}

bool my_print1(const GearEngineContext&);

int main()
{
 double initValues[] = { 0, 1, 0, 0};
 GearUserProc uf1(diffun1), uf2(j1);

 Gear solver;
 solver.set(2, &uf1, &uf2); // task dimention is 2

 GearAssistant& assistant = solver.assistant;
 assistant.params.setParameter("HINI", 0.0001);  // starting step of integration 
 assistant.params.setParameter("EPS", 0.001);    // step error constant 

 GearEngineCallbackEmbed printHandler(my_print);
 solver.setCallback(printHandler);

 try
 {
  Message rc=solver.solve(0, 5, initValues);     // integrate for t from 0 to 5
  printf("\n\nsteps: %d, derCalls: %d, jacCalls: %d\n\n",
    assistant.statistics.steps, assistant.statistics.derCalls,
    assistant.statistics.jacCalls);
 }
 catch(Message& e)
 {
   printf("\n\n%s",e.getMessageID().data());
 }
 printf("\n");
}

bool my_print1(const GearEngineContext& context)
{
 printf("\n%10.5f:",context.t);
 // print y values and their derivatives
 int n = context.n;
 for(int i=0; i < n; ++i) printf("\t%10.5f %10.5f,",
     context.Z[i], context.Z[i+n]/context.h);
 printf("\t(%10.5f  %d)\n", context.h, context.q);
 return true;
}