NAG Library Manual, Mark 27.2
Interfaces:  FL   CL   CPP   AD
/* D02BJ_P0W_F C++ Header Example Program.
*
* Copyright 2021 Numerical Algorithms Group.
*
* Mark 27.2, 2021.
*
*/

#include <iostream>
#include <math.h>
#include <nag.h>
#include <nagad.h>
#include <stdio.h>
using namespace std;

#ifdef __cplusplus
extern "C"
{
#endif
static void NAG_CALL fcn(void *&       ad_handle,
const double &x,
const double  y[],
double        f[],
Integer       iuser[],
double        ruser[]);
static void NAG_CALL g(void *&       ad_handle,
const double &x,
const double  y[],
double &      retval,
Integer       iuser[],
double        ruser[]);
#ifdef __cplusplus
}
#endif

int main(void)
{
int           exit_status = 0;
const Integer n           = 3;

cout << "D02BJ_P0W_F C++ Header Example Program Results\n\n";

Integer iw, iuser[1];
double  pi, tol, x, xinit, xend, ruser[2];
double *y = 0, *w = 0;

iw = 20 * n;
w  = new double[iw];
y  = new double[n];

xinit = 0.0;
xend  = 10.0;
pi    = nag_math_pi;

cout << "\n\nCase: no intermediate output, root-finding\n";
tol = 1.0e-5;
cout << "\n  Calculation with tol = " << tol << endl;

x    = xinit;
y[0] = 0.5;
y[1] = 0.5;
y[2] = pi / 5.0;

const double alpha = -0.032;
const double beta  = -0.02;
ruser[0]           = alpha;
ruser[1]           = beta;

Integer ifail     = 0;
void *  ad_handle = 0;
nag::ad::d02bj(ad_handle, x, xend, n, y, fcn, tol, "D", nullptr, g, w,
-1, iuser, -1, ruser, ifail);

cout.setf(ios::fixed);
cout.setf(ios::right);
cout.precision(3);
cout << "\n  Root of Y(1) = 0.0 at ";
cout.width(5);
cout << x << endl;
cout << "\n  Solution is  ";
cout.precision(4);
for (int i = 0; i < 3; ++i)
{
cout.width(10);
cout << y[i];
}
cout << endl;

delete[] w;
delete[] y;
return exit_status;
}

static void NAG_CALL fcn(void *&       ad_handle,
const double &x,
const double  y[],
double        f[],
Integer       iuser[],
double        ruser[])
{
double alpha, beta;
alpha = ruser[0];
beta  = ruser[1];

f[0] = tan(y[2]);
f[1] = alpha * tan(y[2]) / y[1] + beta * y[1] / cos(y[2]);
f[2] = alpha / (y[1] * y[1]);
}

static void NAG_CALL g(void *&       ad_handle,
const double &x,
const double  y[],
double &      retval,
Integer       iuser[],
double        ruser[])
{
retval = y[0];
}