NAG Library Manual, Mark 30
```/* nag::ad::d02pr Tangent over Tangent Example Program.
*/

#include <dco.hpp>
#include <iostream>

// Function which calls NAG AD routines.
template <typename T> void func(T &eps, std::vector<T> &y);

// Driver with the tangent calls.
// Solves the problem x''=x/r^3, y''=-y/r^3, with r=sqrt(x^2+y^2) and
// initial conditions x(0)=1-eps, x'(0)=0, y(0)=0, y'(0)=sqrt((1+eps)/(1-eps))
// by solving the ODE system:
//    y1'=y3, y2'=y4, y3'=-y1/r^3, y4'=-y2/r^3
// over the range [0,6*pi].
// Also computes the sum of all Hessian elements d2y/deps2.
void driver(const double &epsv, std::vector<double> &yv, double &d2y_deps2);

int main()
{
std::cout
<< "nag::ad::d02pr Tangent over Tangent Example Program Results\n\n";

// Parameter epsilon
double epsv = 0.7;
// Solution y
std::vector<double> yv;
// Sum of Hessian elements
double d2y_deps2;

// Call driver
driver(epsv, yv, d2y_deps2);

// Print outputs
std::cout << "\n Derivatives calculated: Second order tangents\n";
std::cout << " Computational mode    : algorithmic\n";

std::cout.setf(std::ios::scientific, std::ios::floatfield);
std::cout.precision(6);
std::cout << "\n Solution computed with required tolerance " << 1e-4
<< std::endl;
for (std::size_t i = 0; i < yv.size(); i++)
{
std::cout << " y" << i + 1 << " = " << yv[i] << std::endl;
}

// Print derivatives
std::cout
<< "\n Sum of all Hessian elements of solution vector y w.r.t. problem parameter epsilon :\n";
std::cout << " sum_i [d^2 y/deps^2]_i = " << d2y_deps2 << std::endl;

return 0;
}

// Driver with the tangent calls.
// Solves the problem x''=x/r^3, y''=-y/r^3, with r=sqrt(x^2+y^2) and
// initial conditions x(0)=1-eps, x'(0)=0, y(0)=0, y'(0)=sqrt((1+eps)/(1-eps))
// by solving the ODE system:
//    y1'=y3, y2'=y4, y3'=-y1/r^3, y4'=-y2/r^3
// over the range [0,6*pi].
// Also computes the sum of all hessian elements d2y/deps2.
void driver(const double &epsv, std::vector<double> &yv, double &d2y_deps2)
{
using mode = dco::gt1s<dco::gt1s<double>::type>;
using T    = mode::type;

// Variable to differentiate w.r.t.
T eps                            = epsv;
dco::derivative(dco::value(eps)) = 1.0;
dco::value(dco::derivative(eps)) = 1.0;

// Solution y, variable to differentiate
std::vector<T> y;

// Call the NAG AD Lib functions
func(eps, y);

// Extract the computed solutions
yv = dco::passive_value(y);

// Tangents of y
d2y_deps2 = 0;
for (T const &yi : y)
{
d2y_deps2 += dco::derivative(dco::derivative(yi));
}
}

// function which calls NAG AD Library routines
template <typename T> void func(T &eps, std::vector<T> &y)
{
// Active variables
const Integer n = 4, npts = 6;
const Integer liwsav = 130, lrwsav = 350 + 32 * n;

std::vector<T>       thresh(n, 1e-10), ypnow(n), rwsav(lrwsav);
std::vector<Integer> iwsav(liwsav);

// Set parameters for the integrator.
Integer method = -3;
T tol = 1e-4, hstart = 0.0, tnow, tend = 6.0 * nag_math_pi, tstart = 0.0;
T tinc  = (tend - tstart) / ((double)npts);
T twant = tstart + tinc;
// Set initial conditions
y.resize(n);
y[0] = 1.0 - eps;
y[1] = y[2] = 0.0;
y[3]        = sqrt((1.0 + eps) / (1.0 - eps));
// Create AD configuration data object
Integer           ifail = 0;

// Initialize Runge-Kutta method for integrating ODE
ifail = 0;
method, hstart, iwsav.data(), rwsav.data(), ifail);

const T &          t,
const Integer &    n,
const T            y[],
T                  yp[])
{
T r   = 1.0 / sqrt(y[0] * y[0] + y[1] * y[1]);
r     = r * r * r;
yp[0] = y[2];
yp[1] = y[3];
yp[2] = -y[0] * r;
yp[3] = -y[1] * r;
};
do
{
ifail = 0;
// Solve an initial value problem for a 1st order system of ODEs