NAG Library Manual, Mark 28.5
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.5, 2022.
*
*/

#include <math.h>
#include <nag.h>
#include <stdio.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL out(Integer neq, double *xsol, const double y[],
Nag_User *comm);
static void NAG_CALL fcn(Integer neq, double x, const double y[], double f[],
Nag_User *comm);
static double NAG_CALL g(Integer neq, double x, const double y[],
Nag_User *comm);
#ifdef __cplusplus
}
#endif

struct user {
double xend, h;
Integer k;
Integer *use_comm;
};

int main(void) {
static Integer use_comm[2] = {1, 1};
Integer exit_status = 0, i, j, neq;
Nag_User comm;
double pi, tol, x, y[3];
struct user s;
NagError fail;

INIT_FAIL(fail);

/* For communication with user-supplied functions
* assign address of user defined structure
* to Nag pointer.
*/
s.use_comm = use_comm;
comm.p = (Pointer)&s;

neq = 3;
s.xend = 10.0;
/* nag_math_pi (x01aac).
* pi
*/
pi = nag_math_pi;
printf("\nCase 1: intermediate output, root-finding\n");
for (j = 4; j <= 5; ++j) {
tol = pow(10.0, (double)(-j));
printf("\n  Calculation with tol = %10.1e\n", tol);
x = 0.0;
y[0] = 0.5;
y[1] = 0.5;
y[2] = pi / 5.0;
s.k = 4;
s.h = (s.xend - x) / (double)(s.k + 1);
printf("\n     X         Y(1)         Y(2)         Y(3)\n");

* Ordinary differential equation solver using a
* variable-order variable-step Adams method (Black Box)
*/
nag_ode_ivp_adams_zero_simple(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, out,
g, &comm, &fail);
if (fail.code != NE_NOERROR) {
fail.message);
exit_status = 1;
goto END;
}

printf("\n  Root of Y(1) = 0.0 at %7.3f\n", x);
printf("\n  Solution is");
for (i = 0; i < 3; ++i)
printf("%10.5f", y[i]);
printf("\n");
}
printf("\n\nCase 2: no intermediate output, root-finding\n");
for (j = 4; j <= 5; ++j) {
tol = pow(10.0, (double)(-j));
printf("\n  Calculation with tol = %10.1e\n", tol);
x = 0.0;
y[0] = 0.5;
y[1] = 0.5;
y[2] = pi / 5.0;

/* nag_ode_ivp_adams_zero_simple (d02cjc), see above. */
nag_ode_ivp_adams_zero_simple(neq, fcn, &x, y, s.xend, tol, Nag_Mixed,
NULLFN, g, &comm, &fail);
if (fail.code != NE_NOERROR) {
fail.message);
exit_status = 1;
goto END;
}
printf("\n  Root of Y(1) = 0.0 at %7.3f\n", x);
printf("\n  Solution is");
for (i = 0; i < 3; ++i)
printf("%10.5f", y[i]);
printf("\n");
}
printf("\n\nCase 3: intermediate output, no root-finding\n");
for (j = 4; j <= 5; ++j) {
tol = pow(10.0, (double)(-j));
printf("\n  Calculation with tol = %10.1e\n", tol);
x = 0.0;
y[0] = 0.5;
y[1] = 0.5;
y[2] = pi / 5.0;
s.k = 4;
s.h = (s.xend - x) / (double)(s.k + 1);
printf("\n     X         Y(1)         Y(2)         Y(3)\n");

/* nag_ode_ivp_adams_zero_simple (d02cjc), see above. */
nag_ode_ivp_adams_zero_simple(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, out,
NULLDFN, &comm, &fail);
if (fail.code != NE_NOERROR) {
fail.message);
exit_status = 1;
goto END;
}
}

printf("\n\nCase 4: no intermediate output, no root-finding");
printf(" ( integrate to xend)\n");
for (j = 4; j <= 5; ++j) {
tol = pow(10.0, (double)(-j));
printf("\n  Calculation with tol = %10.1e\n", tol);
x = 0.0;
y[0] = 0.5;
y[1] = 0.5;
y[2] = pi / 5.0;
printf("\n     X         Y(1)         Y(2)         Y(3)\n");
printf("%8.2f", x);
for (i = 0; i < 3; ++i)
printf("%13.5f", y[i]);
printf("\n");

/* nag_ode_ivp_adams_zero_simple (d02cjc), see above. */
nag_ode_ivp_adams_zero_simple(neq, fcn, &x, y, s.xend, tol, Nag_Mixed,
NULLFN, NULLDFN, &comm, &fail);
if (fail.code != NE_NOERROR) {
fail.message);
exit_status = 1;
goto END;
}

printf("%8.2f", x);
for (i = 0; i < 3; ++i)
printf("%13.5f", y[i]);
printf("\n");
}
END:
return exit_status;
}

static void NAG_CALL out(Integer neq, double *xsol, const double y[],
Nag_User *comm) {
Integer i;
struct user *s = (struct user *)comm->p;

printf("%8.2f", *xsol);
for (i = 0; i < 3; ++i) {
printf("%13.5f", y[i]);
}
printf("\n");
*xsol = s->xend - (double)s->k * s->h;
s->k--;
}

static void NAG_CALL fcn(Integer neq, double x, const double y[], double f[],
Nag_User *comm) {
double pwr;
struct user *s = (struct user *)comm->p;

if (s->use_comm[0]) {
printf("(User-supplied callback fcn, first invocation.)\n");
s->use_comm[0] = 0;
}

f[0] = tan(y[2]);
f[1] = -0.032 * tan(y[2]) / y[1] - 0.02 * y[1] / cos(y[2]);

pwr = y[1];
f[2] = -0.032 / (pwr * pwr);
}

static double NAG_CALL g(Integer neq, double x, const double y[],
Nag_User *comm) {
struct user *s = (struct user *)comm->p;

if (s->use_comm[1]) {
printf("(User-supplied callback g, first invocation.)\n");
s->use_comm[1] = 0;
}

return y[0];
}