// sim_difeq.java  solve simultaneous differential equations
//                by Runge-Kutta by Gills method
//
// solve system of equations Y_i = Y'_i(X, Y_1, ..., Y_n) for i=1,n
//
// user sets N number of equations (in N dependent variables)
//           Y array of N dependent variables with initial values
//           F array of N first derivatives of Y's (computed in loop)
//           Q array of N for use by Runge
//           X independent variable with initial value
//           XLIM limit value of X with value
//           H step size, delta X with value
//           M setp counter, initialize to zero
//           K returned 1=compute derivatives, 2=have Y values at X+H 

public strictfp class sim_difeq
{
  static int M;
  static double X;

  public static void  main(String args[])
  {
    int N = 4;                   // number of equations
    double Y[] = new double[N];  // dependent variables
    double F[] = new double[N];  // computed derivatives, Y'
    double Q[] = new double[N];  // temporaries for Runge
    double XLIM, H;
    int K;

    System.out.println("sim_difeq.java solving "+N+" equations");
    X = 0.0;      // initial independent variable
    XLIM = 1.0;   // final value of independent variable
    H = 0.1;      // step size
    M = 0;        // step counter
                  // initial condition for dependent variables
    for(int i=0; i<N; i++) Y[i] = 0.0; 
	

    System.out.println("simple numeric check, answer 1, 2, 3, 4");
    while(X < XLIM) // simple numeric check
    {
      K = 1;
      while(K == 1)
      {
	K = Runge(N, Y, F, Q, H); // X, M global, K returned
        F[0] = 1.0;
        F[1] = 2.0;
        F[2] = 3.0;
        F[3] = 4.0;
      } // X incremented by H
      System.out.println("X="+X+" Y[0]="+Y[0]+" Y[1]="+Y[1]);
      System.out.println("                    Y[2]="+Y[2]+" Y[3]="+Y[3]);
    } // X at XLIM (done)

    System.out.println(" ");
    System.out.println("fourth order check 1, 1/2, 1/6, 1/24");
    X = 0.0;      // initial independent variable
    XLIM = 1.0;   // final value of independent variable
    H = 0.01;     // step size
    M = 0;        // step counter
                  // initial condition for dependent variables
    for(int i=0; i<N; i++) Y[i] = 0.0; 
	
    while(X < XLIM) // fourth order check
    {
      K = 1;
      while(K == 1)
      {
        K = Runge(N, Y, F, Q, H);
	F[0] = 1.0;
	F[1] = Y[0];
	F[2] = Y[1];
	F[3] = Y[2];
      } // X incremented by H
    } // X at XLIM (done)
    System.out.println("X="+X+" Y[0]="+Y[0]+"  Y[1]="+Y[1]);
    System.out.println("                    Y[2]="+Y[2]+" Y[3]="+Y[3]);
    System.out.println("Y[3]*24.0 ="+Y[3]*24.0);
    System.out.println(" ");
    System.out.println("sim_difeq.java finished");
  } // end main of class sim_difeq

  static int Runge(int N, double Y[], double F[], double Q[],
                   double H) // Gills method
  {
    double A;
    int K;

    M++;     // integration step
    K = 1;   // exit condition
    if(M == 1)
    {
	for(int i=0; i<N; i++) Q[i] = 0.0;
    }
    else if(M == 2)
    {
      A = 0.5;
      X = X+0.5*H;
      for(int i=0; i<N; i++)
      {
        Y[i] = Y[i] + A*F[i]*H;
        Q[i] = 2.0*A*H*F[i];
      }
    }
    else if(M == 3)
    {
      A = 0.2928932188134524756;
      for(int i=0; i<N; i++)
      {
	Y[i] = Y[i] + A*(F[i]*H-Q[i]);
        Q[i] = 2.0*A*H*F[i] + (1.0-3.0*A)*Q[i];
      }
    }
    else if(M == 4)
    {
      A = 1.7071067811865475244;
      X = X+0.5*H;
      for(int i=0; i<N; i++)
      {
	Y[i] = Y[i] + A*(F[i]*H-Q[i]);
        Q[i] = 2.0*A*H*F[i] + (1.0-3.0*A)*Q[i];
      }
    }
    else if(M == 5)
    {
      A = 0.2928932188134524756;
      for(int i=0; i<N; i++)
      {
        Y[i] = Y[i] + H*F[i]/6.0 - Q[i]/3.0;
      }
      M = 0;
      K = 2;
    }
    else
    {
      System.out.println("ERROR in Runge, bad M="+M);
      M = 0;
      K = 2;
    }
    return K;
  } // end Runge
} // end sim_difeq.java