/* nag_ode_bvp_ps_lin_cgl_deriv (d02udc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.1, 2023.
 */
#include <math.h>
#include <nag.h>

#ifdef __cplusplus
extern "C" {
#endif
static double NAG_CALL fcn(double x);
static double NAG_CALL deriv(double x);
#ifdef __cplusplus
}
#endif

int main(void) {
  /*  Scalars */
  Integer exit_status = 0;
  Integer i, n;
  double a = 0.0, b = 1.5, scale;
  double teneps = 100.0 * nag_machine_precision;
  double uxerr = 0.0;
  /*  Arrays */
  double *f = 0, *fd = 0, *x = 0;
  /* NAG types */
  Nag_Boolean reqerr = Nag_FALSE;
  NagError fail;

  INIT_FAIL(fail);
  printf("nag_ode_bvp_ps_lin_cgl_deriv (d02udc) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%*[^\n] ", &n);
  if (!(f = NAG_ALLOC((n + 1), double)) || !(fd = NAG_ALLOC((n + 1), double)) ||
      !(x = NAG_ALLOC((n + 1), double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* nag_ode_bvp_ps_lin_cgl_grid (d02ucc).
   * Generate Chebyshev Gauss-Lobatto solution grid.
   */
  nag_ode_bvp_ps_lin_cgl_grid(n, a, b, x, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_ps_lin_cgl_grid (d02ucc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Evaluate the function on Chebyshev grid. */
  for (i = 0; i < n + 1; i++)
    f[i] = fcn(x[i]);

  /* nag_ode_bvp_ps_lin_cgl_deriv (d02udc).
   * Differentiate a function using function values on Chebyshev grid.
   */
  nag_ode_bvp_ps_lin_cgl_deriv(n, f, fd, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_ps_lin_cgl_deriv (d02udc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  scale = 2.0 / (b - a);
  for (i = 0; i < n + 1; i++)
    fd[i] = scale * fd[i];

  /* Print function and its derivative. */
  printf("Original function f and numerical derivative fx\n\n");
  printf("%8s%11s%11s\n", "x", "f", "fx");
  for (i = 0; i < n + 1; i++)
    printf("%10.4f %10.4f %10.4f\n", x[i], f[i], fd[i]);

  if (reqerr) {
    for (i = 0; i < n + 1; i++)
      uxerr = MAX(uxerr, fabs(fd[i] - deriv(x[i])));
    printf("fx is within a multiple %8" NAG_IFMT " of machine precision.\n",
           100 * ((Integer)(uxerr / teneps) + 1));
  }
END:
  NAG_FREE(f);
  NAG_FREE(fd);
  NAG_FREE(x);
  return exit_status;
}

static double NAG_CALL fcn(double x) { return 2.0 * x + exp(-x); }

static double NAG_CALL deriv(double x) { return 2.0 - exp(-x); }