/* nag_ode_ivp_rk_range (d02pcc) 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>
#include <nagx01.h>

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

#define NEQ 2
#define ZERO 0.0
#define ONE 1.0
#define TWO 2.0
#define FOUR 4.0

int main(void)
{
  Integer exit_status=0, i, j, neq, nout;
  NagError fail;
  Nag_ErrorAssess errass;
  Nag_ODE_RK opt;
  Nag_RK_method method;
  Nag_User comm;
  double hstart, pi, tend, tgot, *thres=0, tinc, tol, tstart, twant, *ygot=0;
  double *ymax=0, *ypgot=0, *ystart=0;

  INIT_FAIL(fail);

  Vprintf("nag_ode_ivp_rk_range (d02pcc) Example Program Results\n");

  /* Set initial conditions and input for nag_ode_ivp_rk_setup (d02pvc) */
  neq = NEQ;
  if (neq>=1)
    {
      if ( !( thres = NAG_ALLOC(neq, double)) ||
           !( ygot = NAG_ALLOC(neq, double)) ||
           !( ymax = NAG_ALLOC(neq, double)) ||
           !( ypgot = NAG_ALLOC(neq, double)) ||
           !( ystart = NAG_ALLOC(neq, double)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      exit_status = 1;
      return exit_status;
    }

  /* nag_pi (x01aac).
   * pi
   */
  pi = X01AAC;
  tstart = ZERO;
  ystart[0] = ZERO;
  ystart[1] = ONE;
  tend = TWO*pi;
  for (i=0; i<neq; i++)
    thres[i] = 1.0e-8;
  errass = Nag_ErrorAssess_off;
  hstart = ZERO;
  method  = Nag_RK_2_3;

  /*
   * Set control for output
   */
  nout = 8;
  tinc = (tend-tstart)/nout;

  for (i=1; i<=2; i++)
    {
      if (i==1)
        tol = 1.0e-3;
      else
        tol = 1.0e-4;
      /* nag_ode_ivp_rk_setup (d02pvc).
       * Setup function for use with nag_ode_ivp_rk_range (d02pcc)
       * and/or nag_ode_ivp_rk_onestep (d02pdc)
       */
      nag_ode_ivp_rk_setup(neq, tstart, ystart, tend, tol, thres, method,
                           Nag_RK_range, errass, hstart, &opt, &fail);
      if (fail.code != NE_NOERROR)
        {
          Vprintf("Error from nag_ode_ivp_rk_setup (d02pvc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }
      Vprintf("\nCalculation with tol = %8.1e\n\n",tol);
      Vprintf("    t          y1         y2\n\n");
      Vprintf("%8.3f   %8.3f   %8.3f\n", tstart, ystart[0], ystart[1]);
      for (j=nout-1; j>=0; j--)
        {
          twant = tend - j*tinc;
          /* nag_ode_ivp_rk_range (d02pcc).
           * Ordinary differential equations solver, initial value
           * problems over a range using Runge-Kutta methods
           */
          nag_ode_ivp_rk_range(neq, f, twant, &tgot, ygot, ypgot, ymax, &opt,
                               &comm, &fail);
          if (fail.code != NE_NOERROR)
            {
              Vprintf("Error from nag_ode_ivp_rk_range (d02pcc).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }

          Vprintf("%8.3f   %8.3f   %8.3f\n", tgot, ygot[0], ygot[1]);
        }
      Vprintf("\nCost of the integration in evaluations of f is"
              " %ld\n\n", opt.totfcn);
      /* nag_ode_ivp_rk_free (d02ppc).
       * Freeing function for use with the Runge-Kutta suite (d02p
       * functions)
       */
      nag_ode_ivp_rk_free(&opt);
    }
 END:
  if (thres) NAG_FREE(thres);
  if (ygot) NAG_FREE(ygot);
  if (ymax) NAG_FREE(ymax);
  if (ypgot) NAG_FREE(ypgot);
  if (ystart) NAG_FREE(ystart);
  return exit_status;
}
static void f(Integer neq, double t, double y[], double yp[], Nag_User *comm)

{
  yp[0] = y[1];
  yp[1] = -y[0];
}