/*
 *  Eular3tst.c: Modified Eular Method  修正オイラー法
 */
#include <stdio.h>
#include <math.h>

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

int Eular3(double (*f)(), double *t, double *x, double *x0, double h)
/*
    Differencial equation
    by using modified Eular method.
    Result is stored in "x".
 */
{
    double x1;

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

/*  Main  */
int main()
{
    double x, x0, h, t, a, error, tmax;
    int i, 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; x0 = x/(1-h);
/*  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);
        }
    /*  -------------------------- */
        Eular3(f, &t, &x, &x0, h);
    /*  -------------------------- */
    }
    return 0;
}

/*  ---------   実行例: 教科書 p.149 例2   ---------  */
/*
Z:\src\Pascal>Eular3tst
 0.000  1.0000000e+00  0.0000000e+00
 0.250  7.8030579e-01  1.5050064e-03
 0.500  6.0593777e-01 -5.9288645e-04
 0.750  4.7430115e-01  1.9346019e-03
 1.000  3.6641095e-01 -1.4684940e-03
 1.250  2.8925996e-01  2.7551653e-03
 1.500  2.2033177e-01 -2.7983865e-03
 1.750  1.7799021e-01  4.2162697e-03
 2.000  1.3044494e-01 -4.8903471e-03
 2.250  1.1211422e-01  6.7149914e-03
 2.500  7.3826136e-02 -8.2588627e-03
 2.750  7.4831279e-02  1.0903418e-02
 3.000  3.6035592e-02 -1.3751477e-02
 3.250  5.6632945e-02  1.7858737e-02
 3.500  7.4353882e-03 -2.2761995e-02
 3.750  5.2877999e-02  2.9360254e-02
 4.000 -1.9268530e-02 -3.7584169e-02
 4.250  6.2608800e-02  4.8344566e-02
 4.500 -5.0886600e-02 -6.1995597e-02
 4.750  8.8307042e-02  7.9655347e-02
 5.000 -9.5482540e-02 -1.0222049e-01
 5.250  1.3652668e-01  1.3127916e-01
 5.500 -1.6442988e-01 -1.6851665e-01
 5.750  2.1956541e-01  2.1638263e-01
 6.000 -2.7531261e-01 -2.7779136e-01

*/