#include <graphics.h>
#include <conio.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

float matkyv(float T, float X, float xPrime)
{
  return -.001*sin(X);
}

void InitialCond2ndOrder(float LowerLimit,
			 float UpperLimit,
			 float InitialValue,
			 float InitialDeriv,
			   int NumIntervals,
			 float Function(float,float,float))

/*----------------------------------------------------------------------------*/
/*-                                                                          -*/
/*-    Input:  LowerLimit, UpperLimit, InitialValue, InitialDeriv,           -*/
/*-            NumIntervals, NumReturn                                       -*/
/*-                                                                          -*/
/*-  Purpose:  This unit solves a 2nd order ordinary differential            -*/
/*-            equation with specified initial conditions.                   -*/
/*-            Given a function of the form                                  -*/
/*-                       X" = TNTargetF(T, X, X')                           -*/
/*-            and initial conditions                                        -*/
/*-                 X[LowerLimit] = InitialValue,                            -*/
/*-                X'[LowerLimit] = InitialDeriv                             -*/
/*-            the Runge-Kutta method of two variables (X and X') is used to -*/
/*-            solve for X and X' in the interval [LowerLimit, UpperLimit].  -*/
/*-            The interval is divided up into NumIntervals sub-intervals.   -*/
/*-                                                                          -*/
/*-  User-defined Functions: TNTargetF(T, X, xPrime: real) : real;           -*/
/*-                                                                          -*/
/*-                                                                          -*/
/*-    Global Variables:  LowerLimit   : real;    Lower limit of integration -*/
/*-                       UpperLimit   : real;    Upper limit of integration -*/
/*-                       InitialValue : real;    Value of X at LowerLimit   -*/
/*-                       InitialDeriv : real;    Deriv. of X at LowerLimit  -*/
/*-                                               wrong                      -*/
/*-                                                                          -*/
/*-               Notes:  NumIntervals should be an integer                  -*/
/*-                       multiple of NumReturn if you want                  -*/
/*-                       evenly spaced values returned.                     -*/
/*-                                                                          -*/
/*-                                                                          -*/
/*----------------------------------------------------------------------------*/


{
  float Spacing, HalfSpacing;           /* Size of each subinterval  */
  int Index;
  float F1x, F2x, F3x, F4x,
  F1xPrime, F2xPrime,
  F3xPrime, F4xPrime;                   /* Iteration variables  */
  float T, X, xPrime;                   /* Values pushed and pulled  */
					/* from the stack.           */

  /* procedure InitialCond2ndOrder */
  T = LowerLimit;
  X = InitialValue;
  xPrime = InitialDeriv;
//  fprintf(output,"\n%e  %e  %e\n",T,X,xPrime);
  printf("%e  %e  %e\n",T,X,xPrime);
  Spacing = (UpperLimit-LowerLimit)/NumIntervals;
  HalfSpacing = Spacing / 2;
  for (Index = 1;Index <= NumIntervals;Index++)
  {
    F1x = Spacing * xPrime;
    F1xPrime = Spacing * Function(T, X, xPrime);
    F2x = Spacing * (xPrime + 0.5 * F1xPrime);
    F2xPrime = Spacing * Function(T + HalfSpacing, X + 0.5 * F1x,
		xPrime + 0.5 * F1xPrime);
    F3x = Spacing * (xPrime + 0.5 * F2xPrime);
    F3xPrime = Spacing * Function(T + HalfSpacing, X + 0.5 * F2x,
		xPrime + 0.5 * F2xPrime);
    F4x = Spacing * (xPrime + F3xPrime);
    F4xPrime = Spacing * Function(T + Spacing, X + F3x, xPrime +
					F3xPrime);
    X = X + (F1x + 2 * F2x + 2 * F3x + F4x) / 6;
    xPrime = xPrime + (F1xPrime + 2 * F2xPrime + 2 * F3xPrime + F4xPrime) / 6;
    T = T + Spacing;
//    fprintf(output,"%e  %e  %e\n",T,X,xPrime);
    printf("%e  %e  %e\n",T,X,xPrime);
//    lineto(T/1.5,X*100+240);
  }
} /* procedure InitialCond2ndOrder */



int main()
{
     /* request auto detection */
    int gdriver = DETECT, gmode, errorcode;

    /* initialize graphics mode */
    initgraph(&gdriver, &gmode, "C:\\BORLANDC\\BGI");

    /* read result of initialization */
    errorcode = graphresult();

    if (errorcode != grOk)  /* an error occurred */
    {
       printf("Graphics error: %s\n", grapherrormsg(errorcode));
       printf("Press any key to halt:");
       getch();
       exit(1);             /* return with error code */
    }
    moveto(0,240+2*100);
    setcolor(RED);
    InitialCond2ndOrder(0.0,900.0,2.0,0.0,350,&matkyv);
    getch();
    closegraph();
    return 0;
}