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

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

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

#define NEQ 2
#define MNP 70

int main(void)
{

  double a, b, c[NEQ][NEQ], d[NEQ][NEQ];
  Integer i, j;
  double eps;
  double  x[MNP], y[NEQ][MNP];
  Integer neq, mnp, np;
  double gam[NEQ], tol;
  Nag_User comm;

  Vprintf("d02gbc Example Program Results\n");

  /* For communication with function fcnf()
   * assign address of eps to comm.p.
   */
  comm.p = (Pointer)&eps;

  neq = NEQ;
  mnp = MNP;
  tol = 1.0e-3;
  np = 0;
  a = 0.0;
  b = 1.0;

  for (i=0; i<neq; ++i)
    {
      gam[i] = 0.0;
      for (j=0; j<neq; ++j)
        {
          c[i][j] = 0.0;
          d[i][j] = 0.0;
        }
    }
  c[0][0] = 1.0;
  d[1][0] = 1.0;
  gam[1] = 1.0;
  for (i=1; i<=2; ++i)
    {
      eps = pow (10.0, (double) -i);
      Vprintf("\nProblem with epsilon = %7.4f\n", eps);
      d02gbc(neq, fcnf, NULLFN, a, b, &c[0][0], &d[0][0], gam,
             mnp, &np, x, &y[0][0], tol, &comm, NAGERR_DEFAULT);

      Vprintf("\nApproximate solution on final mesh of %ld points\n", np);
      Vprintf("    X(I)        Y(1,I)\n");
      for (j=0; j<np; ++j)
        Vprintf("%9.4f   %9.4f\n", x[j], y[0][j]);

    }
  return EXIT_SUCCESS;
}


static void fcnf(Integer neq, double x, double f[], Nag_User *comm)
{
#define F(I, J) f[(I)*neq+J]
  double *eps = (double *)comm->p;

  F(0,0) = 0.0;
  F(0,1) = 1.0;
  F(1,0) = 0.0;
  F(1,1) = -1.0/ *eps;

}