/* 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 <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


#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);

  Vprintf("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)) )
        {
          Vprintf("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)
    {
      Vprintf("\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)
            {
              Vprintf("Error from nag_ode_bvp_fd_nonlin_fixedbc"
                      " (d02gac).\n%s\n", fail.message);
              exit_status = 2;
            }
          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;
        }
      else
        {
          Vprintf("Error from nag_ode_bvp_fd_nonlin_fixedbc (d02gac).\n%s\n",
                  fail.message);
          exit_status = 1;
        }
    }
 END:
  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 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]);
}