/*************************************************************************/
/* Copyright (C) 1995-1997 Marlis Hochbruck, Christian Lubich, Hubert    */
/* Selhofer. All rights reserved.                                        */
/*                                                                       */
/* C code written by Mathias Froehlich.                                  */
/*                                                                       */
/* This code is part of a copyrighted package. For details, see the file */
/* "COPYING" in the top-level directory.                                 */
/*************************************************************************/

/*****************************************************************/
/* The Van Der Pol's Equation.                                   */
/* 2 - dimensional Pooblem.                                      */
/* for details see:                                              */
/* Hairer, Noersett, Wanner:                                     */
/* Solving Ordinary Differential Equations II, p 157.            */
/* Springer-Verlag 1991.                                         */
/*****************************************************************/
#include <math.h>
#include "ftypes.h"
#include "exp4.h"

#define T0        0.0
#define TEND      2.0
#define NTSPAN    201

#define NEQ       2
#define EPS       1e-6
#define MU        (1.0/EPS)

void vdpex(finteger neq,fdouble *fy,fdouble *y)
{
  fy[0] = y[1];
  fy[1] = MU*((1.0-y[0]*y[0])*y[1]-y[0]);
}

void vdpjac(finteger neq,fdouble *jac,finteger ld,fdouble *y)
{
#define Jac(I,J) jac[(I) +(J)*ld]
  Jac(0,0)=0.0;
  Jac(0,1)=1.0;
  Jac(1,0)= MU - MU*2.0*y[0]*y[1];
  Jac(1,1)= MU*(1.0-y[0]*y[0]);
#undef Jac
}

int main()
{
  int i;
  fdouble t=T0;
  fdouble y0[NEQ];
  FILE *file;
  fdouble Tspan[NTSPAN];
  finteger NTspan = NTSPAN;
  fdouble  yout[NEQ*NTSPAN];
  finteger ldyout = NEQ;

  printf("\n###############################################################\n"
	"Exp4Ex Example: The Van'd Pol equation, 2 - dimensional Problem\n"
	"###############################################################\n"
	"\nThis example contains some output features!\n\n"
	 "Press return please ...\n");
  getchar();

  /* set initial values */
  y0[0]=2.0;
  y0[1]=0.0;

  /* Documentation of the functions below in ../exp4/exp4.h */

  /* initializes Tspan to write output every 1e2 seconds */
  for (i=0;i<NTSPAN;Tspan[i++]=(fdouble)(i*(T0-TEND)/(NTSPAN-1)));

  /* opening an output file */
  file = fopen("vdpex.out","w");

  /* output into the file 'vdpex.out' at Tspan points */
  DExp4SetOutput(DenseFile,file,Tspan,NTspan);

  /* output into file 'vdpex.out' at all computed points */
/*   DExp4SetOutput(AllFile,file); */

  /* output into file 'vdpex.out' only at the end point */
/*   DExp4SetOutput(EndFile,file); */

  /* dense output into the array yout at tspan points.
     (here 'vdpex.out' stays empty) */
/*   DExp4SetOutput(Dense,yout,ldyout,Tspan,NTspan); */

  /* call exp4ex */
  DExp4Ex(&t,TEND,y0,NEQ,&vdpex,&vdpjac,Exp4False);

  /* comment out the following if you commented out the
     'Dense output line above' */
/*   printf("The result of dense output into the output array\n" */
/* 	"Press return please ...\n\n"); */
/*   getchar(); */
/*   for (i=0;i<NTSPAN;i++) */
/*     printf("t = %e, y = [%e %e]\n",Tspan[i],yout[2*i],yout[2*i+1]); */

  fclose(file);

  return 1;
}