/* Runge-Kutta numerical integration
 * of a system of ordinary differential equations.
 *
 * Reference:
 * Thomas, Benku, "The Runge-Kutta Methods," BYTE 11, #4, April 1986
 * This program is adapted from code supplied with that article.
 */

/*include "mconf.h"*/
#define UNK 1
#include "prec.h"
#include "int.h"
#define ORDER 7
#define ERRTERM 0

#if !(ORDER-4)
/*
Fourth order Runge-Kutta method specified by
Felberg, E., COMPUTING,	6(1970)p61-71
*/
#define INUM 6
#define HMAXL ((DOUBLE )0.5)
#define HLIM1 ((DOUBLE )3.0)

#if UNK
DOUBLE al[5] = {
 2.50000000000000000000E-1, /* 1.0/4.0 */
 3.75000000000000000000E-1, /* 3.0/8.0 */
 9.23076923076923076923E-1, /* 12.0/13.0 */
 1.00000000000000000000E0, /* 1.0 */
 5.00000000000000000000E-1, /* 1.0/2.0 */
};

DOUBLE b[15] = {
 2.50000000000000000000E-1, /* 1.0/4.0 */ /* b(1,1) */

 9.37500000000000000000E-2, /* 3.0/32.0 */ /* b(1,2) */
 2.81250000000000000000E-1, /* 9.0/32.0 */ /* b(2,2) */

 8.79380974055530268548E-1, /* 1932.0/2197.0 */ /* b(1,3) */
-3.27719617660446062813E0, /* -7200.0/2197.0 */
 3.32089212562585343650E0, /* 7296.0/2197.0 */ /* b(3,3) */

 2.03240740740740740741E0, /* 439.0/216.0 */ /* b(1,4) */
-8.00000000000000000000E0, /* -8.0 */
 7.17348927875243664717E0, /* 3680.0/513.0 */
-2.05896686159844054581E-1, /* -845.0/4104.0 */

-2.96296296296296296296E-1, /* -8.0/27.0 */ /* b(1,5) */
 2.00000000000000000000E0, /* 2.0 */
-1.38167641325536062378E0, /* -3544.0/2565.0 */
 4.52972709551656920078E-1, /* 1859.0/4104.0 */
-2.75000000000000000000E-1, /* -11.0/40.0 */ /* b(5,5) */
};

#if ERRTERM
DOUBLE er[6] = {
 2.77777777777777777778E-3, /* 1.0/360.0 */
 0.00000000000000000000E0, /* 0.0 */
-2.99415204678362573099E-2, /* -384.0/12825.0 */
-2.91998936735778841042E-2, /* -41743.0/1429560.0 */
 2.00000000000000000000E-2, /* 1.0/50.0 */
 3.63636363636363636364E-2, /* 2.0/55.0 */
};
#endif

DOUBLE a[6] = {
 1.18518518518518518519E-1, /* 16.0/135.0 */
 0.00000000000000000000E0, /* 0.0 */
 5.18986354775828460039E-1, /* 6656.0/12825.0 */
 5.06131490342016657806E-1, /* 28561.0/56430.0 */
-1.80000000000000000000E-1, /* -9.0/50.0 */
 3.63636363636363636364E-2, /* 2.0/55.0 */
};
#endif /* UNK */

#if DEC
static short al[20] = {
0037600,0000000,0000000,0000000,
0037700,0000000,0000000,0000000,
0040154,0047304,0166116,0142354,
0040200,0000000,0000000,0000000,
0040000,0000000,0000000,0000000,
};
static short b[60] = {
0037600,0000000,0000000,0000000,
0037300,0000000,0000000,0000000,
0037620,0000000,0000000,0000000,
0040141,0017434,0106113,0005137,
0140521,0136625,0004105,0164734,
0040524,0104577,0020106,0154200,
0040402,0011366,0102275,0120457,
0141000,0000000,0000000,0000000,
0040745,0106471,0061516,0054324,
0137522,0153224,0132645,0026551,
0137627,0132045,0166411,0075502,
0040400,0000000,0000000,0000000,
0140260,0155305,0150113,0006654,
0037747,0166011,0175402,0077301,
0137614,0146314,0146314,0146315,
};
#if ERRTERM
static short er[24] = {
0036066,0005540,0133013,0060266,
0000000,0000000,0000000,0000000,
0136765,0043753,0063165,0135345,
0136757,0032235,0105775,0077355,
0036643,0153412,0036560,0121727,
0037024,0171011,0047440,0112362,
};
#endif
static short a[24] = {
0037362,0134726,0044017,0025635,
0000000,0000000,0000000,0000000,
0040004,0156112,0026325,0015221,
0040001,0110725,0053176,0122300,
0137470,0050753,0102436,0134122,
0037024,0171011,0047440,0112362,
};
#endif /* DEC */

#if IBMPC
static short al[20] = {
0x0000,0x0000,0x0000,0x3fd0,
0x0000,0x0000,0x0000,0x3fd8,
0xd89e,0x9d89,0x89d8,0x3fed,
0x0000,0x0000,0x0000,0x3ff0,
0x0000,0x0000,0x0000,0x3fe0,
};
static short b[60] = {
0x0000,0x0000,0x0000,0x3fd0,
0x0000,0x0000,0x0000,0x3fb8,
0x0000,0x0000,0x0000,0x3fd2,
0x614c,0x9189,0x23e3,0x3fec,
0xbd3c,0xa108,0x37b2,0xc00a,
0xdb10,0xe408,0x912f,0x400a,
0xb426,0xd097,0x425e,0x4000,
0x0000,0x0000,0x0000,0xc020,
0xcb1a,0x2c69,0xb1a7,0x401c,
0xa5ad,0x96b4,0x5ad2,0xbfca,
0x2f68,0xbda1,0xf684,0xbfd2,
0x0000,0x0000,0x0000,0x4000,
0x61b6,0xba09,0x1b58,0xbff6,
0x4fd8,0x3f60,0xfd81,0x3fdc,
0x999a,0x9999,0x9999,0xbfd1,
};
#if ERRTERM
static short er[24] = {
0x6c17,0x16c1,0xc16c,0x3f66,
0x0000,0x0000,0x0000,0x0000,
0xb75d,0x6cce,0xa8fd,0xbf9e,
0xafde,0xb17f,0xe693,0xbf9d,
0x147b,0x47ae,0x7ae1,0x3f94,
0x129e,0x29e4,0x9e41,0x3fa2,
};
#endif
static short a[24] = {
0xe574,0xc901,0x573a,0x3fbe,
0x0000,0x0000,0x0000,0x0000,
0xa352,0x459a,0x9b89,0x3fe0,
0xd498,0xaacf,0x323a,0x3fe0,
0xd70a,0x70a3,0x0a3d,0xbfc7,
0x129e,0x29e4,0x9e41,0x3fa2,
};
#endif /* IBMPC */

#if MIEEE
static short al[20] = {
0x3fd0,0x0000,0x0000,0x0000,
0x3fd8,0x0000,0x0000,0x0000,
0x3fed,0x89d8,0x9d89,0xd89e,
0x3ff0,0x0000,0x0000,0x0000,
0x3fe0,0x0000,0x0000,0x0000,
};
static short b[60] = {
0x3fd0,0x0000,0x0000,0x0000,
0x3fb8,0x0000,0x0000,0x0000,
0x3fd2,0x0000,0x0000,0x0000,
0x3fec,0x23e3,0x9189,0x614c,
0xc00a,0x37b2,0xa108,0xbd3c,
0x400a,0x912f,0xe408,0xdb10,
0x4000,0x425e,0xd097,0xb426,
0xc020,0x0000,0x0000,0x0000,
0x401c,0xb1a7,0x2c69,0xcb1a,
0xbfca,0x5ad2,0x96b4,0xa5ad,
0xbfd2,0xf684,0xbda1,0x2f68,
0x4000,0x0000,0x0000,0x0000,
0xbff6,0x1b58,0xba09,0x61b6,
0x3fdc,0xfd81,0x3f60,0x4fd8,
0xbfd1,0x9999,0x9999,0x999a,
};
#if ERRTERM
static short er[24] = {
0x3f66,0xc16c,0x16c1,0x6c17,
0x0000,0x0000,0x0000,0x0000,
0xbf9e,0xa8fd,0x6cce,0xb75d,
0xbf9d,0xe693,0xb17f,0xafde,
0x3f94,0x7ae1,0x47ae,0x147b,
0x3fa2,0x9e41,0x29e4,0x129e,
};
#endif
static short a[24] = {
0x3fbe,0x573a,0xc901,0xe574,
0x0000,0x0000,0x0000,0x0000,
0x3fe0,0x9b89,0x459a,0xa352,
0x3fe0,0x323a,0xaacf,0xd498,
0xbfc7,0x0a3d,0x70a3,0xd70a,
0x3fa2,0x9e41,0x29e4,0x129e,
};
#endif /* MIEEE */
#endif

#if !(ORDER-5)
/*
Fifth order Runge-Kutta method specified by
Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 5)
*/
#define INUM 8
#define HMAXL ((DOUBLE )1.0)
#define HLIM1 ((DOUBLE )4.0)

#if UNK
DOUBLE al[7] = {
 5.55555555555555555556E-2, /* 1.0/18.0 */
 1.66666666666666666667E-1, /* 1.0/6.0 */
 2.22222222222222222222E-1, /* 2.0/9.0 */
 6.66666666666666666667E-1, /* 2.0/3.0 */
 1.00000000000000000000E0, /* 1.0 */
 8.88888888888888888889E-1, /* 8.0/9.0 */
 1.00000000000000000000E0, /* 1.0 */
};

DOUBLE b[28] = {
 5.55555555555555555556E-2, /* 1.0/18.0 */ /* b(1,1) */

-8.33333333333333333333E-2, /* -1.0/12.0 */ /* b(1,2) */
 2.50000000000000000000E-1, /* 1.0/4.0 */

-2.46913580246913580247E-2, /* -2.0/81.0 */ /* b(1,3) */
 1.48148148148148148148E-1, /* 4.0/27.0 */
 9.87654320987654320988E-2, /* 8.0/81.0 */

 1.21212121212121212121E0, /* 40.0/33.0 */ /* b(1,4) */
-3.63636363636363636364E-1, /* -4.0/11.0 */
-5.09090909090909090909E0, /* -56.0/11.0 */
 4.90909090909090909091E0, /* 54.0/11.0 */

-5.05479452054794520548E0, /* -369.0/73.0 */ /* b(1,5) */
 9.86301369863013698630E-1, /* 72.0/73.0 */
 2.45662100456621004566E1, /* 5380.0/219.0 */
-2.10359589041095890411E1, /* -12285.0/584.0 */
 1.53824200913242009132E0, /* 2695.0/1752.0 */

-9.78226711560044893378E0, /* -8716.0/891.0 */ /* b(1,6) */
 2.20875420875420875421E0, /* 656.0/297.0 */
 4.43546576879910213244E1, /* 39520.0/891.0 */
-3.78181818181818181818E1, /* -416.0/11.0 */
 1.92592592592592592593E0, /* 52.0/27.0 */
 0.00000000000000000000E0, /* 0.0 */

 1.17773437500000000000E1, /* 3015.0/256.0 */ /* b(1,7) */
-2.25000000000000000000E0, /* -9.0/4.0 */
-5.40897435897435897436E1, /* -4219.0/78.0 */
 4.67578125000000000000E1, /* 5985.0/128.0 */
-1.40364583333333333333E0, /* -539.0/384.0 */
 0.00000000000000000000E0, /* 0.0 */
 2.08233173076923076923E-1, /* 693.0/3328.0 */
};

#if ERRTERM
DOUBLE er[8] = {
 5.15625000000000000000E-2, /* 33.0/640.0 */
 0.00000000000000000000E0, /* 0.0 */
-4.06153846153846153846E-1, /* -132.0/325.0 */
 3.97767857142857142857E-1, /* 891.0/2240.0 */
-1.03125000000000000000E-1, /* -33.0/320.0 */
-1.04285714285714285714E-1, /* -73.0/700.0 */
 1.07091346153846153846E-1, /* 891.0/8320.0 */
 5.71428571428571428571E-2, /* 2.0/35.0 */
};
#endif

DOUBLE a[8] = {
 8.90625000000000000000E-2, /* 57.0/640.0 */
 0.00000000000000000000E0, /* 0.0 */
-2.46153846153846153846E-1, /* -16.0/65.0 */
 6.14732142857142857143E-1, /* 1377.0/2240.0 */
 3.78125000000000000000E-1, /* 121.0/320.0 */
 0.00000000000000000000E0, /* 0.0 */
 1.07091346153846153846E-1, /* 891.0/8320.0 */
 5.71428571428571428571E-2, /* 2.0/35.0 */
};
#endif /* UNK */

#if DEC
static short al[28] = {
0037143,0107070,0161616,0034344,
0037452,0125252,0125252,0125253,
0037543,0107070,0161616,0034344,
0040052,0125252,0125252,0125253,
0040200,0000000,0000000,0000000,
0040143,0107070,0161616,0034344,
0040200,0000000,0000000,0000000,
};

static short b[112] = {
0037143,0107070,0161616,0034344,

0137252,0125252,0125252,0125253,
0037600,0000000,0000000,0000000,

0136712,0042607,0163267,0047403,
0037427,0132045,0166411,0075502,
0037312,0042607,0163267,0047403,

0040233,0023311,0131154,0115447,
0137672,0027213,0121350,0135057,
0140642,0164272,0027213,0121351,
0040635,0013505,0150564,0056427,

0140641,0140340,0070070,0016016,
0040174,0077077,0017617,0143744,
0041304,0103631,0020746,0044172,
0141250,0044644,0151151,0032232,
0040304,0162435,0034507,0047122,

0141034,0102052,0102743,0153707,
0040415,0056072,0116340,0015625,
0041461,0065453,0061214,0144437,
0141427,0042721,0072135,0013506,
0040366,0102275,0120457,0064114,
0000000,0000000,0000000,0000000,

0041074,0070000,0000000,0000000,
0140420,0000000,0000000,0000000,
0141530,0055745,0137133,0162676,
0041473,0004000,0000000,0000000,
0140263,0125252,0125252,0125253,
0000000,0000000,0000000,0000000,
0037525,0035423,0130473,0011661,
};

#if ERRTERM
static short er[8] = {
0037123,0031463,0031463,0031463,
0000000,0000000,0000000,0000000,
0137717,0171545,0116300,0031152,
0037713,0124072,0101650,0035204,
0137323,0031463,0031463,0031463,
0137325,0111677,0121140,0106157,
0037333,0051265,0025522,0132453,
0037152,0007240,0165016,0120352,
};
#endif

static short a[32] = {
0037266,0063146,0063146,0063146,
0000000,0000000,0000000,0000000,
0137574,0007700,0176017,0140374,
0040035,0057425,0170537,0012761,
0037701,0114631,0114631,0114632,
0000000,0000000,0000000,0000000,
0037333,0051265,0025522,0132453,
0037152,0007240,0165016,0120352,
};
#endif /* DEC */

#if IBMPC
static short al[28] = {
0xc71c,0x1c71,0x71c7,0x3fac,
0x5555,0x5555,0x5555,0x3fc5,
0xc71c,0x1c71,0x71c7,0x3fcc,
0x5555,0x5555,0x5555,0x3fe5,
0x0000,0x0000,0x0000,0x3ff0,
0xc71c,0x1c71,0x71c7,0x3fec,
0x0000,0x0000,0x0000,0x3ff0,
};

static short b[112] = {
0xc71c,0x1c71,0x71c7,0x3fac,

0x5555,0x5555,0x5555,0xbfb5,
0x0000,0x0000,0x0000,0x3fd0,

0xe9e0,0xfcd6,0x48b0,0xbf99,
0x2f68,0xbda1,0xf684,0x3fc2,
0xe9e0,0xfcd6,0x48b0,0x3fb9,

0x9365,0x364d,0x64d9,0x3ff3,
0x1746,0x745d,0x45d1,0xbfd7,
0x745d,0x45d1,0x5d17,0xc014,
0x8ba3,0xba2e,0xa2e8,0x4013,

0x0382,0x0e07,0x381c,0xc014,
0xf8fc,0xe3f1,0x8fc7,0x3fef,
0xc90f,0x243c,0x90f3,0x4038,
0x2693,0x9a4d,0x0934,0xc035,
0xe9ca,0xa728,0x9ca3,0x3ff8,

0x7af9,0x50bc,0x9085,0xc023,
0x0373,0x539c,0xab87,0x4001,
0x9924,0x6c51,0x2d65,0x4046,
0xa2e9,0x2e8b,0xe8ba,0xc042,
0xed09,0xb425,0xd097,0x3ffe,
0x0000,0x0000,0x0000,0x0000,

0x0000,0x0000,0x8e00,0x4027,
0x0000,0x0000,0x0000,0xc002,
0x7cb8,0xb7cb,0x0b7c,0xc04b,
0x0000,0x0000,0x6100,0x4047,
0x5555,0x5555,0x7555,0xbff6,
0x0000,0x0000,0x0000,0x0000,
0x6276,0x7627,0xa762,0x3fca,
};

#if ERRTERM
static short er[8] = {
0x6666,0x6666,0x6666,0x3faa,
0x0000,0x0000,0x0000,0x0000,
0x064d,0xb398,0xfe6c,0xbfd9,
0x0750,0x5075,0x7507,0x3fd9,
0x6666,0x6666,0x6666,0xbfba,
0x118e,0xf44c,0xb277,0xbfba,
0x56a5,0xa56a,0x6a56,0x3fbb,
0xd41d,0x1d41,0x41d4,0x3fad,
};
#endif

static short a[32] = {
0xcccd,0xcccc,0xcccc,0x3fb6,
0x0000,0x0000,0x0000,0x0000,
0xf820,0x1f81,0x81f8,0xbfcf,
0xe2be,0xbe2b,0xabe2,0x3fe3,
0x3333,0x3333,0x3333,0x3fd8,
0x0000,0x0000,0x0000,0x0000,
0x56a5,0xa56a,0x6a56,0x3fbb,
0xd41d,0x1d41,0x41d4,0x3fad,
};
#endif /* IBMPC */


#if MIEEE
static short al[28] = {
0x3fac,0x71c7,0x1c71,0xc71c,
0x3fc5,0x5555,0x5555,0x5555,
0x3fcc,0x71c7,0x1c71,0xc71c,
0x3fe5,0x5555,0x5555,0x5555,
0x3ff0,0x0000,0x0000,0x0000,
0x3fec,0x71c7,0x1c71,0xc71c,
0x3ff0,0x0000,0x0000,0x0000,
};

static short b[112] = {
0x3fac,0x71c7,0x1c71,0xc71c,

0xbfb5,0x5555,0x5555,0x5555,
0x3fd0,0x0000,0x0000,0x0000,

0xbf99,0x48b0,0xfcd6,0xe9e0,
0x3fc2,0xf684,0xbda1,0x2f68,
0x3fb9,0x48b0,0xfcd6,0xe9e0,

0x3ff3,0x64d9,0x364d,0x9365,
0xbfd7,0x45d1,0x745d,0x1746,
0xc014,0x5d17,0x45d1,0x745d,
0x4013,0xa2e8,0xba2e,0x8ba3,

0xc014,0x381c,0x0e07,0x0382,
0x3fef,0x8fc7,0xe3f1,0xf8fc,
0x4038,0x90f3,0x243c,0xc90f,
0xc035,0x0934,0x9a4d,0x2693,
0x3ff8,0x9ca3,0xa728,0xe9ca,

0xc023,0x9085,0x50bc,0x7af9,
0x4001,0xab87,0x539c,0x0373,
0x4046,0x2d65,0x6c51,0x9924,
0xc042,0xe8ba,0x2e8b,0xa2e9,
0x3ffe,0xd097,0xb425,0xed09,
0x0000,0x0000,0x0000,0x0000,

0x4027,0x8e00,0x0000,0x0000,
0xc002,0x0000,0x0000,0x0000,
0xc04b,0x0b7c,0xb7cb,0x7cb8,
0x4047,0x6100,0x0000,0x0000,
0xbff6,0x7555,0x5555,0x5555,
0x0000,0x0000,0x0000,0x0000,
0x3fca,0xa762,0x7627,0x6276,
};

#if ERRTERM
static short er[8] = {
0x3faa,0x6666,0x6666,0x6666,
0x0000,0x0000,0x0000,0x0000,
0xbfd9,0xfe6c,0xb398,0x064d,
0x3fd9,0x7507,0x5075,0x0750,
0xbfba,0x6666,0x6666,0x6666,
0xbfba,0xb277,0xf44c,0x118e,
0x3fbb,0x6a56,0xa56a,0x56a5,
0x3fad,0x41d4,0x1d41,0xd41d,
};
#endif

static short a[8] = {
0x3fb6,0xcccc,0xcccc,0xcccd,
0x0000,0x0000,0x0000,0x0000,
0xbfcf,0x81f8,0x1f81,0xf820,
0x3fe3,0xabe2,0xbe2b,0xe2be,
0x3fd8,0x3333,0x3333,0x3333,
0x0000,0x0000,0x0000,0x0000,
0x3fbb,0x6a56,0xa56a,0x56a5,
0x3fad,0x41d4,0x1d41,0xd41d,
};
#endif /* MIEEE */
#endif

#if !(ORDER-7)
/*
Seventh order Runge-Kutta method specified by
Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 7)
*/
#if UNK
DOUBLE al[12] = {
 2.50000000000000000000E-1, /* 1.0/4.0 */
 8.33333333333333333333E-2, /* 1.0/12.0 */
 1.25000000000000000000E-1, /* 1.0/8.0 */
 4.00000000000000000000E-1, /* 2.0/5.0 */
 5.00000000000000000000E-1, /* 1.0/2.0 */
 8.57142857142857142857E-1, /* 6.0/7.0 */
 1.42857142857142857143E-1, /* 1.0/7.0 */
 6.66666666666666666667E-1, /* 2.0/3.0 */
 2.85714285714285714286E-1, /* 2.0/7.0 */
 1.00000000000000000000E0, /* 1.0 */
 3.33333333333333333333E-1, /* 1.0/3.0 */
 1.00000000000000000000E0, /* 1.0 */
};

DOUBLE b[78] = {
 2.50000000000000000000E-1, /* 1.0/4.0 */ /* b(1,1) */

 6.94444444444444444444E-2, /* 5.0/72.0 */ /* b(1,2) */
 1.38888888888888888889E-2, /* 1.0/72.0 */

 3.12500000000000000000E-2, /* 1.0/32.0 */ /* b(1,3) */
 0.00000000000000000000E0, /* 0.0 */
 9.37500000000000000000E-2, /* 3.0/32.0 */

 8.48000000000000000000E-1, /* 106.0/125.0 */ /* b(1,4) */
 0.00000000000000000000E0, /* 0.0 */
-3.26400000000000000000E0, /* -408.0/125.0 */
 2.81600000000000000000E0, /* 352.0/125.0 */

 2.08333333333333333333E-2, /* 1.0/48.0 */ /* b(1,5) */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 2.42424242424242424242E-1, /* 8.0/33.0 */
 2.36742424242424242424E-1, /* 125.0/528.0 */

-5.26030820491461890879E-1, /* -1263.0/2401.0 */ /* b(1,6) */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 1.51209723221385028965E0, /* 39936.0/26411.0 */
-2.42796562038544545833E0, /* -64125.0/26411.0 */
 2.29904206580591420242E0, /* 5520.0/2401.0 */

 9.43877551020408163265E-2, /* 37.0/392.0 */ /* b(1,7) */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 1.72725340136054421769E-1, /* 1625.0/9408.0 */
-1.33333333333333333333E-1, /* -2.0/15.0 */
 9.07738095238095238095E-3, /* 61.0/6720.0 */

 6.73172643543013913384E-1, /* 17176.0/25515.0 */ /* b(1,8) */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
-1.84612972761120909269E0, /* -47104.0/25515.0 */
 2.62896825396825396825E0, /* 1325.0/504.0 */
-1.63793846756809719773E0, /* -41792.0/25515.0 */
 1.38799725651577503429E-1, /* 20237.0/145800.0 */
 7.09794238683127572016E-1, /* 4312.0/6075.0 */

-1.32355962793280577537E-1, /* -23834.0/180075.0 */ /* b(1,9) */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
-3.92886802216248280388E-2, /* -77824.0/1980825.0 */
-1.00437160021708126664E0, /* -636635.0/633864.0 */
 8.46473969179508538109E-1, /* 254048.0/300125.0 */
-2.61428571428571428571E-2, /* -183.0/7000.0 */
 7.27272727272727272727E-1, /* 8.0/11.0 */
-8.58733103631062814736E-2, /* -324.0/3773.0 */

 1.67539473684210526316E0, /* 12733.0/7600.0 */ /* b(1,10) */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
-3.83387559808612440191E0, /* -20032.0/5225.0 */
 5.68786134370015948963E0, /* 456485.0/80256.0 */
-5.97880701754385964912E0, /* -42599.0/7125.0 */
 3.71959429824561403509E-1, /* 339227.0/912000.0 */
-2.46172248803827751196E-1, /* -1029.0/4180.0 */
 1.20809659090909090909E0, /* 1701.0/1408.0 */
 2.11554276315789473684E0, /* 5145.0/2432.0 */

-1.32573976092494611013E-1, /* -27061.0/204120.0 */ /* b(1,11) */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 1.44114870040795966722E-1, /* 40448.0/280665.0 */
-1.13049726765004542782E0, /* -1353775.0/1197504.0 */
 6.92220262590632961003E-1, /* 17662.0/25515.0 */
-6.14600480109739368999E-2, /* -71687.0/1166400.0 */
 4.35555555555555555556E-1, /* 98.0/225.0 */
 6.25000000000000000000E-2, /* 1.0/16.0 */
 3.23473936899862825789E-1, /* 3773.0/11664.0 */
 0.00000000000000000000E0, /* 0.0 */

 1.29066820276497695853E0, /* 11203.0/8680.0 */ /* b(1,12) */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
-3.19597821533305404273E0, /* -38144.0/11935.0 */
 5.13725605711492808267E0, /* 2354425.0/458304.0 */
-5.16411674347158218126E0, /* -84046.0/16275.0 */
 4.11356915933528836755E-1, /* 673309.0/1636800.0 */
 5.51788856304985337243E-1, /* 4704.0/8525.0 */
 8.68493401759530791789E-1, /* 9477.0/10912.0 */
-1.03729838709677419355E0, /* -1029.0/992.0 */
 0.00000000000000000000E0, /* 0.0 */
 2.13782991202346041056E0, /* 729.0/341.0 */ /* b(12,12) */
};
#if ERRTERM
DOUBLE er[13] = {
-2.08333333333333333333E-3, /* -1.0/480.0 */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
-4.26666666666666666667E-2, /* -16.0/375.0 */
-4.54734848484848484848E-3, /* -2401.0/528000.0 */
 1.81893939393939393939E-2, /* 2401.0/132000.0 */
 1.72585227272727272727E-2, /* 243.0/14080.0 */
-1.25052083333333333333E-1, /* -2401.0/19200.0 */
-4.22222222222222222222E-2, /* -19.0/450.0 */
 1.38068181818181818182E-1, /* 243.0/1760.0 */
 4.30555555555555555556E-2, /* 31.0/720.0 */
};
#endif
DOUBLE a[13] = {
 4.30555555555555555556E-2, /* 31.0/720.0 */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 2.13333333333333333333E-1, /* 16.0/75.0 */
 2.12209595959595959596E-1, /* 16807.0/79200.0 */
 2.12209595959595959596E-1, /* 16807.0/79200.0 */
 1.38068181818181818182E-1, /* 243.0/1760.0 */
 0.00000000000000000000E0, /* 0.0 */
 0.00000000000000000000E0, /* 0.0 */
 1.38068181818181818182E-1, /* 243.0/1760.0 */
 4.30555555555555555556E-2, /* 31.0/720.0 */
};
#endif /* UNK */

#if DEC
static short al[48] = {
0037600,0000000,0000000,0000000,
0037252,0125252,0125252,0125253,
0037400,0000000,0000000,0000000,
0037714,0146314,0146314,0146315,
0040000,0000000,0000000,0000000,
0040133,0066666,0155555,0133333,
0037422,0044444,0111111,0022222,
0040052,0125252,0125252,0125253,
0037622,0044444,0111111,0022222,
0040200,0000000,0000000,0000000,
0037652,0125252,0125252,0125253,
0040200,0000000,0000000,0000000,
};

static short b[312] = {
0037600,0000000,0000000,0000000,

0037216,0034343,0107070,0161616,
0036543,0107070,0161616,0034344,

0037000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0037300,0000000,0000000,0000000,

0040131,0013207,0025402,0006112,
0000000,0000000,0000000,0000000,
0140520,0162540,0040611,0033514,
0040464,0034530,0010142,0046723,

0036652,0125252,0125252,0125253,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0037570,0037017,0101740,0174076,
0037562,0066233,0023311,0131155,

0140006,0124764,0131262,0164601,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0040301,0106146,0170135,0102164,
0140433,0061711,0164727,0105255,
0040423,0021601,0052460,0170212,

0037301,0047136,0005162,0170124,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0037460,0157351,0056114,0120070,
0137410,0104210,0104210,0104211,
0036424,0134513,0112271,0045625,

0040054,0052412,0154224,0040223,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0140354,0046772,0115041,0135544,
0040450,0040404,0010101,0002020,
0140321,0123767,0135610,0056336,
0037416,0020603,0164570,0016341,
0040065,0132423,0041010,0101444,

0137407,0104122,0047202,0146151,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0137040,0166452,0145172,0062047,
0140200,0107477,0121773,0052640,
0040130,0131204,0117212,0111266,
0136726,0024613,0107173,0073424,
0040072,0027213,0121350,0135057,
0137257,0157130,0116336,0151063,

0040326,0071525,0130520,0046634,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0140565,0057067,0140655,0113415,
0040666,0001365,0145353,0177136,
0140677,0051143,0014056,0043422,
0037676,0070567,0062422,0011235,
0137574,0012223,0173456,0046354,
0040232,0121350,0135056,0105643,
0040407,0062415,0074503,0057121,

0137407,0140570,0167170,0143762,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0037423,0111331,0033305,0020413,
0140220,0132042,0066143,0016144,
0040061,0032530,0156564,0107230,
0137173,0136610,0001573,0145404,
0037737,0000443,0042547,0104654,
0037200,0000000,0000000,0000000,
0037645,0117140,0034057,0141062,
0000000,0000000,0000000,0000000,

0040245,0032235,0116156,0066020,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0140514,0105350,0033145,0057324,
0040644,0062146,0150217,0110633,
0140645,0040161,0140675,0173612,
0037722,0116537,0124725,0175235,
0040015,0041010,0152040,0106502,
0040136,0052625,0062531,0053126,
0140204,0143061,0106143,0014306,
0000000,0000000,0000000,0000000,
0040410,0151064,0106443,0044322,
};
#if ERRTERM
static short er[52] = {
0136010,0104210,0104210,0104211,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0137056,0141476,0017547,0012452,
0136225,0000754,0101521,0175550,
0036625,0000754,0101521,0175550,
0036615,0060676,0153033,0166542,
0137400,0006647,0040332,0072016,
0137054,0170465,0074676,0001107,
0037415,0060676,0153033,0166542,
0037060,0055405,0130133,0002660,
};
#endif
static short a[52] = {
0037060,0055405,0130133,0002660,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0037532,0072015,0123500,0155164,
0037531,0046570,0165054,0161770,
0037531,0046570,0165054,0161770,
0037415,0060676,0153033,0166542,
0000000,0000000,0000000,0000000,
0000000,0000000,0000000,0000000,
0037415,0060676,0153033,0166542,
0037060,0055405,0130133,0002660,
};
#endif /* DEC */


#if IBMPC
static short al[48] = {
0x0000,0x0000,0x0000,0x3fd0,
0x5555,0x5555,0x5555,0x3fb5,
0x0000,0x0000,0x0000,0x3fc0,
0x999a,0x9999,0x9999,0x3fd9,
0x0000,0x0000,0x0000,0x3fe0,
0xb6db,0xdb6d,0x6db6,0x3feb,
0x2492,0x9249,0x4924,0x3fc2,
0x5555,0x5555,0x5555,0x3fe5,
0x2492,0x9249,0x4924,0x3fd2,
0x0000,0x0000,0x0000,0x3ff0,
0x5555,0x5555,0x5555,0x3fd5,
0x0000,0x0000,0x0000,0x3ff0,
};

static short b[312] = {
0x0000,0x0000,0x0000,0x3fd0,

0x1c72,0x71c7,0xc71c,0x3fb1,
0xc71c,0x1c71,0x71c7,0x3f8c,

0x0000,0x0000,0x0000,0x3fa0,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x3fb8,

0x4189,0xe560,0x22d0,0x3feb,
0x0000,0x0000,0x0000,0x0000,
0x26e9,0x0831,0x1cac,0xc00a,
0x49ba,0x020c,0x872b,0x4006,

0x5555,0x5555,0x5555,0x3f95,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x1f08,0xf07c,0x07c1,0x3fcf,
0x364e,0x64d9,0x4d93,0x3fce,

0x5d30,0x9656,0xd53e,0xbfe0,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0xb08f,0xde0b,0x318c,0x3ff8,
0xf156,0x3d3a,0x6c79,0xc003,
0x1e11,0x2aa6,0x6470,0x4002,

0x5e0a,0xc14e,0x29cb,0x3fb8,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x9407,0x2b89,0x1bdd,0x3fc6,
0x1111,0x1111,0x1111,0xbfc1,
0x2973,0x7297,0x9729,0x3f82,

0x8812,0x5b12,0x8aa1,0x3fe5,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x376d,0x5344,0x89bf,0xbffd,
0x2082,0x8208,0x0820,0x4005,
0x0b9c,0xf771,0x34fe,0xbffa,
0x039c,0x7d2f,0xc430,0x3fc1,
0x1064,0x6841,0xb6a2,0x3fe6,

0x598d,0x49d0,0xf10a,0xbfc0,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x4c85,0x594f,0x1da5,0xbfa4,
0x6ab4,0xf47f,0x11e7,0xbff0,
0x5257,0x93d1,0x1650,0x3feb,
0x6ee2,0x71cf,0xc531,0xbf9a,
0x1746,0x745d,0x45d1,0x3fe7,
0xda46,0x139b,0xfbcb,0xbfb5,

0x09b3,0xb62a,0xce6a,0x3ffa,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0xb2e2,0xf835,0xabc6,0xc00e,
0x7fcc,0xb95d,0xc05e,0x4016,
0xc8e2,0x6305,0xea4c,0xc017,
0x4254,0xeca2,0xce2e,0x3fd7,
0xc99e,0x7ee5,0x8292,0xbfcf,
0xd174,0x1745,0x545d,0x3ff3,
0x6bca,0xaf28,0xeca1,0x4000,

0x18fe,0x1dcf,0xf82f,0xbfc0,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0xa421,0x26d8,0x725b,0x3fc2,
0x638c,0x4d8c,0x1684,0xbff2,
0x91d3,0x1bae,0x26ab,0x3fe6,
0x7961,0x006f,0x77b1,0xbfaf,
0xf135,0x68ac,0xe024,0x3fdb,
0x0000,0x0000,0x0000,0x3fb0,
0xf846,0x0705,0xb3cc,0x3fd4,
0x0000,0x0000,0x0000,0x0000,

0xcd82,0xb38d,0xa693,0x3ff4,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0xabdb,0x06cc,0x915d,0xc009,
0xf233,0xda11,0x8c8c,0x4014,
0xbef1,0x3837,0xa80e,0xc014,
0xbf54,0xf53a,0x53ab,0x3fda,
0x11a8,0x1a84,0xa841,0x3fe1,
0x2acb,0xacab,0xcab2,0x3feb,
0x6319,0x318c,0x98c6,0xbff0,
0x0000,0x0000,0x0000,0x0000,
0x691a,0x91a4,0x1a46,0x4001,
};
#if ERRTERM
static short er[52] = {
0x1111,0x1111,0x1111,0xbf61,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0xe2a5,0xc3ec,0xd867,0xbfa5,
0x3f6d,0x906a,0xa03d,0xbf72,
0x3f6d,0x906a,0xa03d,0x3f92,
0x7dac,0xdac3,0xac37,0x3f91,
0x4e82,0xe81b,0x01b4,0xbfc0,
0xc049,0xaf37,0x9e26,0xbfa5,
0x7dac,0xdac3,0xac37,0x3fc1,
0x60b6,0xb60b,0x0b60,0x3fa6,
};
#endif
static short a[52] = {
0x60b6,0xb60b,0x0b60,0x3fa6,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x1b4f,0xb4e8,0x4e81,0x3fcb,
0x9c7f,0x1d45,0x29af,0x3fcb,
0x9c7f,0x1d45,0x29af,0x3fcb,
0x7dac,0xdac3,0xac37,0x3fc1,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x7dac,0xdac3,0xac37,0x3fc1,
0x60b6,0xb60b,0x0b60,0x3fa6,
};
#endif /* IBMPC */


#if MIEEE
static short al[48] = {
0x3fd0,0x0000,0x0000,0x0000,
0x3fb5,0x5555,0x5555,0x5555,
0x3fc0,0x0000,0x0000,0x0000,
0x3fd9,0x9999,0x9999,0x999a,
0x3fe0,0x0000,0x0000,0x0000,
0x3feb,0x6db6,0xdb6d,0xb6db,
0x3fc2,0x4924,0x9249,0x2492,
0x3fe5,0x5555,0x5555,0x5555,
0x3fd2,0x4924,0x9249,0x2492,
0x3ff0,0x0000,0x0000,0x0000,
0x3fd5,0x5555,0x5555,0x5555,
0x3ff0,0x0000,0x0000,0x0000,
};

static short b[312] = {
0x3fd0,0x0000,0x0000,0x0000,

0x3fb1,0xc71c,0x71c7,0x1c72,
0x3f8c,0x71c7,0x1c71,0xc71c,

0x3fa0,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x3fb8,0x0000,0x0000,0x0000,

0x3feb,0x22d0,0xe560,0x4189,
0x0000,0x0000,0x0000,0x0000,
0xc00a,0x1cac,0x0831,0x26e9,
0x4006,0x872b,0x020c,0x49ba,

0x3f95,0x5555,0x5555,0x5555,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x3fcf,0x07c1,0xf07c,0x1f08,
0x3fce,0x4d93,0x64d9,0x364e,

0xbfe0,0xd53e,0x9656,0x5d30,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x3ff8,0x318c,0xde0b,0xb08f,
0xc003,0x6c79,0x3d3a,0xf156,
0x4002,0x6470,0x2aa6,0x1e11,

0x3fb8,0x29cb,0xc14e,0x5e0a,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x3fc6,0x1bdd,0x2b89,0x9407,
0xbfc1,0x1111,0x1111,0x1111,
0x3f82,0x9729,0x7297,0x2973,

0x3fe5,0x8aa1,0x5b12,0x8812,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0xbffd,0x89bf,0x5344,0x376d,
0x4005,0x0820,0x8208,0x2082,
0xbffa,0x34fe,0xf771,0x0b9c,
0x3fc1,0xc430,0x7d2f,0x039c,
0x3fe6,0xb6a2,0x6841,0x1064,

0xbfc0,0xf10a,0x49d0,0x598d,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0xbfa4,0x1da5,0x594f,0x4c85,
0xbff0,0x11e7,0xf47f,0x6ab4,
0x3feb,0x1650,0x93d1,0x5257,
0xbf9a,0xc531,0x71cf,0x6ee2,
0x3fe7,0x45d1,0x745d,0x1746,
0xbfb5,0xfbcb,0x139b,0xda46,

0x3ffa,0xce6a,0xb62a,0x09b3,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0xc00e,0xabc6,0xf835,0xb2e2,
0x4016,0xc05e,0xb95d,0x7fcc,
0xc017,0xea4c,0x6305,0xc8e2,
0x3fd7,0xce2e,0xeca2,0x4254,
0xbfcf,0x8292,0x7ee5,0xc99e,
0x3ff3,0x545d,0x1745,0xd174,
0x4000,0xeca1,0xaf28,0x6bca,

0xbfc0,0xf82f,0x1dcf,0x18fe,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x3fc2,0x725b,0x26d8,0xa421,
0xbff2,0x1684,0x4d8c,0x638c,
0x3fe6,0x26ab,0x1bae,0x91d3,
0xbfaf,0x77b1,0x006f,0x7961,
0x3fdb,0xe024,0x68ac,0xf135,
0x3fb0,0x0000,0x0000,0x0000,
0x3fd4,0xb3cc,0x0705,0xf846,
0x0000,0x0000,0x0000,0x0000,

0x3ff4,0xa693,0xb38d,0xcd82,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0xc009,0x915d,0x06cc,0xabdb,
0x4014,0x8c8c,0xda11,0xf233,
0xc014,0xa80e,0x3837,0xbef1,
0x3fda,0x53ab,0xf53a,0xbf54,
0x3fe1,0xa841,0x1a84,0x11a8,
0x3feb,0xcab2,0xacab,0x2acb,
0xbff0,0x98c6,0x318c,0x6319,
0x0000,0x0000,0x0000,0x0000,
0x4001,0x1a46,0x91a4,0x691a,
};
#if ERRTERM
static short er[52] = {
0xbf61,0x1111,0x1111,0x1111,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0xbfa5,0xd867,0xc3ec,0xe2a5,
0xbf72,0xa03d,0x906a,0x3f6d,
0x3f92,0xa03d,0x906a,0x3f6d,
0x3f91,0xac37,0xdac3,0x7dac,
0xbfc0,0x01b4,0xe81b,0x4e82,
0xbfa5,0x9e26,0xaf37,0xc049,
0x3fc1,0xac37,0xdac3,0x7dac,
0x3fa6,0x0b60,0xb60b,0x60b6,
};
#endif
static short a[52] = {
0x3fa6,0x0b60,0xb60b,0x60b6,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x3fcb,0x4e81,0xb4e8,0x1b4f,
0x3fcb,0x29af,0x1d45,0x9c7f,
0x3fcb,0x29af,0x1d45,0x9c7f,
0x3fc1,0xac37,0xdac3,0x7dac,
0x0000,0x0000,0x0000,0x0000,
0x0000,0x0000,0x0000,0x0000,
0x3fc1,0xac37,0xdac3,0x7dac,
0x3fa6,0x0b60,0xb60b,0x60b6,
};
#endif /* MIEEE */

#define INUM 13
#define HMAXL ((DOUBLE )2.0)
#define HLIM1 ((DOUBLE )5.0)
#endif

int inum = INUM;

static DOUBLE *w;
static DOUBLE *au;
static DOUBLE *ays;

/* Initialize pointers in work array.
 * neq >= neqn in rungek() below.  If neq > neqn, unused space
 * is left for extra variables that are not actually integrated.
 */
void rkstart( neq, work )
int neq;
DOUBLE work[];
{
int m;

w = &work[0];
m = neq*inum;
au = &work[m];
m += neq;
ays = &work[m];
}


/* Make a step.
 */
rungek(  neqn, x, yold, h, ynew, delta )
int neqn; /* number of first order equations */
DOUBLE x; /* independent variable */
DOUBLE yold[]; /* initial solution vector at x, size neqn */
DOUBLE h; /* step to take in independent variable */
DOUBLE ynew[]; /* output new solution vector at x+h, size neqn */
DOUBLE delta[]; /* the step, size neqn */
/*DOUBLE err[];*/ /* estimated error, size neqn */
{
DOUBLE *pa, *pal, *pb, *per, *pw, *pb0, *pu, *pys;
DOUBLE xs, usum, esum;
int i, j, m;

/*
 * w[]: neqn blocks of size inum
 *      block of size neqn for temp solution vector ys[]
 *      block of size neqn for output function vector u[] of func()
 */
xs = x;

pys = ays;
pb = yold;
for( i=0; i<neqn; i++ )
	*pys++ = *pb++;  /* ys[i] = yold[i] */

func( xs, ays, au );

pw = w;
pu = au;
for( i=0; i<neqn; i++ )
	{
	*pw = *pu++;  /* w[i*inum] = u[i] */
	pw += inum;
	}

pb0 = (DOUBLE *)&b[0];
pal = (DOUBLE *)&al[0];
for( j=1; j<inum; j++ )
	{
	xs = x + h * *pal++;
	pys = ays;
	for( i=0; i<neqn; i++ )
		{
		usum = 0.0;
		pw = w + i*inum;
		pb = pb0;
		for( m=0; m<j; m++ )
			usum += (*pb++) * (*pw++);
		*pys++ = yold[i] + h * usum;
		}
	func( xs, ays, au );
	pw = w + j;
	pu = au;
	for( i=0; i<neqn; i++ )
		{
		*pw = *pu++;  /* w[j+i*inum] = u[i] */
		pw += inum;
		}
	pb0 += j;
	}

/* evaluate YNEW[i] and the error estimates EST[i] */
pw = w;
for( i=0; i<neqn; i++ )
	{
	pa = (DOUBLE *)&a[0];
#if ERRTERM
	per = (DOUBLE *)&er[0];
#endif
	usum = 0.0;
	esum = 0.0;
	for( j=0; j<inum; j++ )
		{
		if( *pa != 0 )
			{
#if ERRTERM
			esum += (*per++) * (*pw);
#endif
			usum += (*pa++) * (*pw++);
			}
		else
			{ /* skip if coefficient = 0 */
			++pa;
#if ERRTERM
			++per;
#endif
			++pw;
			}
		}
	xs = h * usum; /* the step in this coordinate */
	delta[i] = xs;
	ynew[i] = yold[i] + xs;

#if ERRTERM
	err[i] = h * esum;
#endif
	}
}