/*************************************************************************/
/* 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 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()
{
  fdouble t=T0;
  fdouble y0[NEQ];

  printf("\n###############################################################\n"
	"Exp4Ex Example: The Van'd Pol equation, 2 - dimensional Problem\n"
	"###############################################################\n"
	"\nSimple integrator call!\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 */

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

  return 1;
}