/* nag_ode_bvp_fd_nonlin_fixedbc (d02gac) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.1, 2023.
 *
 */

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

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

#define NEQ 3
#define MNP 40

#define U(I, J) u[(I)*tdu + J]
#define Y(I, J) y[(I)*tdy + J]
#define V(I, J) v[(I)*tdv + J]

int main(void) {
  Integer exit_status = 0, i, j, k, mnp, neq, np, tdu, tdv, tdy, *v = 0;
  NagError fail;
  Nag_User comm;
  double a, b, beta, tol, *u = 0, *x = 0, *y = 0;

  INIT_FAIL(fail);

  printf("nag_ode_bvp_fd_nonlin_fixedbc (d02gac) Example Program Results\n");

  /* For communication with function fcn()
   * assign address of beta to comm.p.
   */
  comm.p = (Pointer)&beta;
  neq = NEQ;
  mnp = MNP;
  tol = 0.001;
  np = 26;
  a = 0.0;
  b = 10.0;
  beta = 0.0;
  if (mnp >= 32 && neq >= 2) {
    if (!(u = NAG_ALLOC(neq * 2, double)) || !(x = NAG_ALLOC(mnp, double)) ||
        !(y = NAG_ALLOC(neq * mnp, double)) ||
        !(v = NAG_ALLOC(neq * 2, Integer))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    tdu = 2;
    tdy = mnp;
    tdv = 2;
  } else {
    exit_status = 1;
    return exit_status;
  }
  for (i = 0; i < neq; ++i)
    for (j = 0; j < 2; ++j) {
      U(i, j) = 0.0;
      V(i, j) = 0;
    }

  V(2, 0) = 1;
  V(0, 1) = 1;
  V(2, 1) = 1;
  U(1, 1) = 1.0;
  U(0, 1) = b;
  x[0] = a;
  for (i = 2; i <= np - 1; ++i)
    x[i - 1] = ((double)(np - i) * a + (double)(i - 1) * b) / (double)(np - 1);
  x[np - 1] = b;
  for (k = 1; k <= 2; ++k) {
    printf("\nProblem with beta = %7.4f\n", beta);
    /* nag_ode_bvp_fd_nonlin_fixedbc (d02gac).
     * Ordinary differential equations solver, for simple
     * nonlinear two-point boundary value problems, using a
     * finite difference technique with deferred correction
     */
    nag_ode_bvp_fd_nonlin_fixedbc(neq, fcn, a, b, u, v, mnp, &np, x, y, tol,
                                  &comm, &fail);

    if (fail.code == NE_NOERROR || fail.code == NE_CONV_ROUNDOFF) {
      if (fail.code == NE_CONV_ROUNDOFF) {
        printf("Error from nag_ode_bvp_fd_nonlin_fixedbc (d02gac)"
               ".\n%s\n",
               fail.message);
        exit_status = 2;
      }
      printf("\nSolution on final mesh of %" NAG_IFMT " points\n", np);
      printf("     X          Y(1)        Y(2)        Y(3)\n");
      for (i = 0; i <= np - 1; ++i) {
        printf(" %9.4f ", x[i]);
        for (j = 0; j < neq; ++j)
          printf(" %9.4f  ", Y(j, i));
        printf("\n");
      }
      beta += 0.2;
    } else {
      printf("Error from nag_ode_bvp_fd_nonlin_fixedbc (d02gac).\n%s\n",
             fail.message);
      exit_status = 1;
    }
  }
END:
  NAG_FREE(u);
  NAG_FREE(x);
  NAG_FREE(y);
  NAG_FREE(v);
  return exit_status;
}

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

  f[0] = y[1];
  f[1] = y[2];
  f[2] = -y[0] * y[2] - *beta * (1.0 - y[1] * y[1]);
}