/* nag_ode_bvp_fd_nonlin_fixedbc (d02gac) Example Program.
 *
 * Copyright 1992 Numerical Algorithms Group.
 *
 * Mark 3, 1992.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

#include <nag.h>
#include <nagx04.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, 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(int argc, char *argv[])
{
  FILE     *fpout;
  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);

  /* Check for command-line IO options */
  fpout = nag_example_file_io(argc, argv, "-results", NULL);

  fprintf(fpout,
          "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)))
        {
          fprintf(fpout, "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)
    {
      fprintf(fpout, "\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)
            {
              fprintf(
                fpout, "Error from nag_ode_bvp_fd_nonlin_fixedbc (d02gac)"
                ".\n%s\n", fail.message);
              exit_status = 2;
            }
          fprintf(fpout, "\nSolution on final mesh of %ld points\n", np);
          fprintf(fpout, "     X          Y(1)        Y(2)        Y(3)\n");
          for (i = 0; i <= np-1; ++i)
            {
              fprintf(fpout, " %9.4f ", x[i]);
              for (j = 0; j < neq; ++j)
                fprintf(fpout, " %9.4f  ", Y(j, i));
              fprintf(fpout, "\n");
            }
          beta += 0.2;
        }
      else
        {
          fprintf(fpout,
                  "Error from nag_ode_bvp_fd_nonlin_fixedbc (d02gac).\n%s\n",
                  fail.message);
          exit_status = 1;
        }
    }
 END:
  if (fpout != stdout) fclose(fpout);
  if (u) NAG_FREE(u);
  if (x) NAG_FREE(x);
  if (y) NAG_FREE(y);
  if (v) NAG_FREE(v);
  return exit_status;
}


static void NAG_CALL fcn(Integer neq, double x, 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]);
}