```/* nag_ode_ivp_bdf_gen (d02ejc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 3, 1992.
* Mark 7 revised, 2001.
* Mark 8 revised, 2004.
*
*/

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

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

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

#define NEQ 3
int main(void)
{
static Integer use_comm = {1, 1, 1, 1};
Integer     exit_status = 0, i, j, neq;
NagError    fail;
Nag_User    comm;
double      tol, x, *y = 0;
struct user s;

INIT_FAIL(fail);

printf("nag_ode_ivp_bdf_gen (d02ejc) Example Program Results\n");

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

neq = NEQ;
if (neq >= 1)
{
if (!(y = NAG_ALLOC(neq, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
else
{
exit_status = 1;
return exit_status;
}
s.xend = 10.0;
printf("\nCase 1: calculating Jacobian internally\n");
printf(" intermediate output, root-finding\n\n");

for (j = 3; j <= 4; ++j)
{
tol = pow(10.0, -(double) j);
printf("\n Calculation with tol = %10.1e\n", tol);
x = 0.0;
y = 1.0;
y = 0.0;
y = 0.0;
s.k = 4;
s.h = (s.xend-x) /(double)(s.k+1);
printf("     X         Y(1)         Y(2)         Y(3)\n");
/* nag_ode_ivp_bdf_gen (d02ejc).
* Ordinary differential equations solver, stiff, initial
* value problems using the Backward Differentiation
* Formulae
*/
nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative,
out, g, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("  Root of Y(1)-0.9 at %5.3f\n", x);
printf("  Solution is ");
for (i = 0; i < 3; ++i)
printf("%7.5f ", y[i]);
printf("\n");
}
printf("\nCase 2: calculating Jacobian by pederv\n");
printf(" intermediate output, root-finding\n\n");

for (j = 3; j <= 4; ++j)
{
tol = pow(10.0, -(double) j);
printf("\n Calculation with tol = %10.1e\n", tol);
x = 0.0;
y = 1.0;
y = 0.0;
y = 0.0;
s.k = 4;
s.h = (s.xend-x) /(double)(s.k+1);
printf("     X         Y(1)         Y(2)         Y(3)\n");
/* nag_ode_ivp_bdf_gen (d02ejc), see above. */
nag_ode_ivp_bdf_gen(neq, fcn, pederv, &x, y, s.xend, tol, Nag_Relative,
out, g, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("  Root of Y(1)-0.9 at %5.3f\n", x);
printf("  Solution is ");
for (i = 0; i < 3; ++i)
printf("%7.5f ", y[i]);
printf("\n");
}
printf("\nCase 3: calculating Jacobian internally\n");
printf(" no intermediate output, root-finding\n\n");
for (j = 3; j <= 4; ++j)
{
tol = pow(10.0, -(double) j);
printf("\n Calculation with tol = %10.1e\n", tol);
x = 0.0;
y = 1.0;
y = 0.0;
y = 0.0;

/* nag_ode_ivp_bdf_gen (d02ejc), see above. */
nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative,
NULLFN, g, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

printf("  Root of Y(1)-0.9 at %5.3f\n", x);
printf("  Solution is ");
for (i = 0; i < 3; ++i)
printf("%7.5f ", y[i]);
printf("\n");
}
printf("\nCase 4: calculating Jacobian internally\n");
printf(" intermediate output, no root-finding\n\n");

for (j = 3; j <= 4; ++j)
{
tol = pow(10.0, -(double) j);
printf("\n Calculation with tol = %10.1e\n", tol);
x = 0.0;
y = 1.0;
y = 0.0;
y = 0.0;
s.k = 4;
s.h = (s.xend-x) /(double)(s.k+1);
printf("     X         Y(1)         Y(2)         Y(3)\n");
/* nag_ode_ivp_bdf_gen (d02ejc), see above. */
nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative,
out, NULLDFN, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
printf("%8.2f", x);
for (i = 0; i < 3; ++i)
printf("%13.5f", y[i]);
printf("\n");
}

printf("\nCase 5: calculating Jacobian internally\n");
printf(
" no intermediate output, no root-finding (integrate to xend)\n\n");

for (j = 3; j <= 4; ++j)
{
tol = pow(10.0, -(double) j);
printf("\n Calculation with tol = %10.1e\n", tol);
x = 0.0;
y = 1.0;
y = 0.0;
y = 0.0;
printf("     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_bdf_gen (d02ejc), see above. */
nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative,
NULLFN, NULLDFN, &comm, &fail);
if (fail.code != NE_NOERROR)
{
printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n",
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:
NAG_FREE(y);
return exit_status;
}

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

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

f = y * -0.04 + y * 1e4 * y;
f = y * 0.04 - y * 1e4 * y - y * 3e7 * y;
f = y * 3e7 * y;
}

static void NAG_CALL pederv(Integer neq, double x, const double y[],
double pw[], Nag_User *comm)
{
#define PW(I, J) pw[((I) -1)*neq + (J) -1]
struct user *s = (struct user *) comm->p;

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

PW(1, 1) = -0.04;
PW(1, 2) = y * 1e4;
PW(1, 3) = y * 1e4;
PW(2, 1) = 0.04;
PW(2, 2) = y * -1e4 - y * 6e7;
PW(2, 3) = y * -1e4;
PW(3, 1) = 0.0;
PW(3, 2) = y * 6e7;
PW(3, 3) = 0.0;
}

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)
{
printf("(User-supplied callback g, first invocation.)\n");
s->use_comm = 0;
}

return y-0.9;
}

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

printf("%8.2f", *xsol);
for (j = 0; j < 3; ++j)
printf("%13.5f", y[j]);
printf("\n");

*xsol = s->xend - (double) s->k * s->h;
s->k--;
}
```