/* nag_ode_bvp_fd_lin_gen (d02gbc) Example Program.

 *

 * Copyright 1992 Numerical Algorithms Group.

 *

 * Mark 3, 1992.

 * Mark 7 revised, 2001.

 * Mark 8 revised, 2004.

 *

 */



#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



#define C(I,J) c[(I)*tdc + J]

#define D(I,J) d[(I)*tdd + J]

#define Y(I,J) y[(I)*tdy + J]



int main(void)

{



  Integer exit_status=0, i, j, mnp, neq, np, tdc, tdd, tdy;

  NagError fail;

  Nag_User comm;

  double a, b, *c=0, *d=0, eps, *gam=0, tol, *x=0, *y=0;



  INIT_FAIL(fail);



  Vprintf("nag_ode_bvp_fd_lin_gen (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;

  if (mnp>=32 && neq>=2)

    {

      if ( !( c = NAG_ALLOC(NEQ*NEQ, double)) ||

           !( d = NAG_ALLOC(NEQ*NEQ, double)) ||

           !( gam = NAG_ALLOC(NEQ, double)) ||

           !( x = NAG_ALLOC(MNP, double)) ||

           !( y = NAG_ALLOC(NEQ*MNP, double)) )

        {

          Vprintf("Allocation failure\n");

          exit_status = -1;

          goto END;

        }

      tdc = neq;

      tdd = neq;

      tdy = mnp;

    }

  else

    {

      exit_status = 1;

      return exit_status;

    }



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

      /* nag_ode_bvp_fd_lin_gen (d02gbc).

       * Ordinary differential equations solver, for general

       * linear two-point boundary value problems, using a finite

       * difference technique with deferred correction

       */

      nag_ode_bvp_fd_lin_gen(neq, fcnf, NULLFN, a, b, c, d, gam,

                             mnp, &np, x, y, tol, &comm, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_ode_bvp_fd_lin_gen (d02gbc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }

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



    }

 END:

  if (c) NAG_FREE(c);

  if (d) NAG_FREE(d);

  if (gam) NAG_FREE(gam);

  if (x) NAG_FREE(x);

  if (y) NAG_FREE(y);

  return exit_status;

}

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;



}