#include <stdio.h>

double fikus(double T, double X)

/*---------------------------------------------------------------*/
/*-         This is the first order differential equation       -*/
/*---------------------------------------------------------------*/

{
  return  T * X;
}
/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-    Input:  LowerLimit, UpperLimit, XInitial, NumReturn, NumIntervals     -*/
/*-   Output:  TValues, XValues, Error                                       -*/
/*-                                                                          -*/
/*-  Purpose:  This unit solves an initial value problem (a first order      -*/
/*-            differential equation with initial condition specified).      -*/
/*-            Given a function of the form                                  -*/
/*-                       dx/dt = TNTargetF(T, X)                            -*/
/*-            and an initial condition                                      -*/
/*-                       X[LowerLimit] = const                              -*/
/*-            the Adams-Bashforth/Adams-Moulton four step predictor/        -*/
/*-            corrector method is used to solve for X in  the interval      -*/
/*-            [LowerLimit, UpperLimit].  The interval is divided up into    -*/
/*-            NumIntervals equal sub-intervals.                             -*/
/*-                                                                          -*/
/*-  User-defined Functions: TNTargetF(T, X : real) : real;                  -*/
/*-                                                                          -*/
/*-      User-defined Types: TNvector = array[0..TNArraySize] of real;       -*/
/*-                                                                          -*/
/*-    Global Variables:  LowerLimit   : real;    Lower limit of integration -*/
/*-                       UpperLimit   : real;    Upper limit of integration -*/
/*-                       XInitial     : real;    Value of X at LowerLimit   -*/
/*-                       NumReturn    : integer; Number of values to return -*/
/*-                       NumIntervals : integer; Number of subintervals     -*/
/*-                       TValues      : TNvector; Value of T between limits -*/
/*-                       XValues      : TNvector; Value of X at TValues     -*/
/*-                       Error        : byte;    Flags if something goes    -*/
/*-                                               wrong                      -*/
/*-                                                                          -*/
/*-               Notes:  NumIntervals should be an integer                  -*/
/*-                       multiple of NumReturn if you want                  -*/
/*-                       evenly spaced values returned.                     -*/
/*-                                                                          -*/
/*-              Errors:  0: No errors                                       -*/
/*-                       1: NumReturn <= 0                                  -*/
/*-                       2: NumIntervals < NumReturn                        -*/
/*-                       3: LowerLimit = UpperLimit                         -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/

typedef double (*function)();   /* Pointers to user defined functions */

void GetFirstThreeValues(double Spacing,
			      double *T,
			      double *X,
			      double *Values,
			      double *DerivValues,
			      double FuncPtr(double,double));  // Runge-Kutta method;

void Adams(double LowerLimit,
		double UpperLimit,
		double XInitial,
		int NumReturn,
		int NumIntervals,
		unsigned char *Error,
		double FuncPtr(double,double))
{
  double Spacing;                 /* Size of each subinterval */
  double T;                       /* Value of T */
  double X, NewX;                 /* Iteration variables */
  int Index;                      /* A counter */
  double DerivValues[5];          /* Derivative at last 4 X values */
  double Values[5];

  *Error = 0;
  if ( NumReturn <= 0 )
    *Error = 1;
  if ( NumIntervals < NumReturn )
    *Error = 2;
  if ( LowerLimit == UpperLimit )
    *Error = 3;
  if ( *Error == 0 )
  {
    for (Index = 0; Index <= 5;Index++) DerivValues[Index] = 0;
    DerivValues[1] = FuncPtr(LowerLimit, XInitial);
    Spacing = (UpperLimit - LowerLimit) / NumIntervals;
    T = LowerLimit;
    X = XInitial;
    printf("%e  %e\n",T,X);
    GetFirstThreeValues(Spacing, &T, &X, Values, DerivValues,FuncPtr);
    for (Index = 4; Index <= NumIntervals; Index++)
    {
      /* Apply predictor */
      NewX = X + Spacing / 24 * (55 * DerivValues[4] - 59 * DerivValues[3] +
		  37 * DerivValues[2] - 9 * DerivValues[1]);

      T = T + Spacing;
      DerivValues[1] = DerivValues[2];
      DerivValues[2] = DerivValues[3];
      DerivValues[3] = DerivValues[4];
      DerivValues[4] = FuncPtr(T, NewX);

      /* Apply corrector */
      NewX = X + Spacing / 24 * (9 * DerivValues[4] + 19 * DerivValues[3] -
	      5 * DerivValues[2] + DerivValues[1]);
      X = NewX;
      printf("%e  %e\n",T,X);     // text output !!!
    }
  }
} /* procedure Adams */

void GetFirstThreeValues(double Spacing,
			      double *T,
			      double *X,
			      double *Values,
			      double *DerivValues,
			      double FuncPtr(double,double))

/*---------------------------------------------------------------*/
/*- Input: Spacing, T, X, Values, DerivValues                   -*/
/*- Output: T, X, Values, DerivValues                           -*/
/*-                                                             -*/
/*- This procedure uses the Runge-Kutta one step method to      -*/
/*- solve for the first three (T + i-Spacing, 1<=i<=3) Values.  -*/
/*- The derivative at these points are also calculated          -*/
/*- (DerivValues).                                              -*/
/*---------------------------------------------------------------*/

{
  int Index;
  double HalfSpacing;
  double F1, F2, F3, F4;

  HalfSpacing = Spacing / 2;
  for (Index = 1; Index < 4; Index++)
  {
    F1 = Spacing * FuncPtr(*T, *X);
    F2 = Spacing * FuncPtr(*T + HalfSpacing, *X + F1 / 2);
    F3 = Spacing * FuncPtr(*T + HalfSpacing, *X + F2 / 2);
    F4 = Spacing * FuncPtr(*T + Spacing, *X + F3);
    *X = *X + (F1 + 2 * F2 + 2 * F3 + F4) / 6;
    *T = *T + Spacing;
    printf("%e  %e\n",*T,*X);     // text output !!!
    DerivValues[Index + 1] = FuncPtr(*T, *X);
  }
} /* procedure GetFirstThreeValues */


void main() {
  int error = 0;

  printf("\n\n\n\n\n\n\n\n\n\n\n\n\n");
  Adams(0,1,1,5,5,&error,&fikus);
}
