/*
 *  Eular2tst.c: Advanced Eular Method  改良オイラー法
 */
#include <stdio.h>
#include <math.h>

double f(double x, double t)
/*  dx/dt = f(x,t) = -x */
{   return -x; }

int Eular2(double (*f)(), double *t, double *x, double h)
/*
    Differencial equation
    by using advanced Eular method.
    Result is stored in "x".
 */
{
    double hh, t1, x1;

    hh = h/2;
    x1 = *x+hh*(*f)(*x,*t); t1 = *t+hh;
    *x +=   h *(*f)(x1,t1); *t +=   h;
    return 0;
}

/*  Main  */
int main()
{
    double x, h, t, a, error, tmax;
    int jstep, nstep, nprnt, nskip;

/*  Parameters */
    h = 0.05; tmax = 6;
    nprnt = 24;
    nstep = (int)(tmax/h+0.5); nskip = nstep/nprnt;
/*  Define initial conditions */
    t = 0; x = 1;
/*  Step on */
    for (jstep=0; jstep<=nstep; jstep++) {
        if (jstep == (jstep/nskip)*nskip) {
            if (t < 700) a = exp(-t); else a = 0;
            error = x-a;
            printf("%6.3f %14.7e %14.7e\n",t, x, error);
        }
    /*  --------------------- */
        Eular2(f, &t, &x, h);
    /*  --------------------- */
    }
    return 0;
}

/*  ---------   実行例: 教科書 p.149 例2   ---------  */
/*
Z:\src\Pascal>Eular2tst
 0.000  1.0000000e+00  0.0000000e+00
 0.250  7.7888502e-01  8.4232614e-05
 0.500  6.0666187e-01  1.3120795e-04
 0.750  4.7251984e-01  1.5328557e-04
 1.000  3.6803862e-01  1.5918050e-04
 1.250  2.8665977e-01  1.5497075e-04
 1.500  2.2327500e-01  1.4483745e-04
 1.750  1.7390555e-01  1.3160655e-04
 2.000  1.3545243e-01  1.1714381e-04
 2.250  1.0550187e-01  1.0264120e-04
 2.500  8.2173822e-02  8.8823744e-05
 2.750  6.4003959e-02  7.6097718e-05
 3.000  4.9851725e-02  6.4656183e-05
 3.250  3.8828761e-02  5.4553427e-05
 3.500  3.0243140e-02  4.5756900e-05
 3.750  2.3555929e-02  3.8182968e-05
 4.000  1.8347360e-02  3.1721103e-05
 4.250  1.4290484e-02  2.6249866e-05
 4.500  1.1130644e-02  2.1647141e-05
 4.750  8.6694916e-03  1.7796373e-05
 5.000  6.7525371e-03  1.4590084e-05
 5.250  5.2594500e-03  1.1931552e-05
 5.500  4.0965068e-03  9.7353195e-06
 5.750  3.1907077e-03  7.9269339e-06
 6.000  2.4851944e-03  6.4422640e-06

*/