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

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

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

#define NEQ 3
#define MNP 40

int main(void)
{
  double a, b;
  Integer i, j, k;
  double u[NEQ][2], x[MNP], y[NEQ][MNP];
  Integer v[NEQ][2];
  Integer np;
  double tol;
  static NagError fail;
  Integer neq, mnp;
  double beta;
  Nag_User comm;

  Vprintf("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;
  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)
    {
      Vprintf("\nProblem with beta = %7.4f\n", beta);
      d02gac(neq, fcn, a, b, &u[0][0], &v[0][0], mnp,
             &np, x, &y[0][0], tol, &comm, &fail);

      if (fail.code == NE_NOERROR || fail.code == NE_CONV_ROUNDOFF)
        {
          Vprintf("\nSolution on final mesh of %ld points\n", np);
          Vprintf("     X          Y(1)        Y(2)        Y(3)\n");
          for (i=0; i<=np-1; ++i)
            {
              Vprintf(" %9.4f ", x[i]);
              for (j=0; j<neq; ++j)
                Vprintf(" %9.4f  ", y[j][i]);
              Vprintf("\n");
            }
          beta += 0.2;
        }
    }
  return EXIT_SUCCESS;
}


static void 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]);
}