/*
 *  Eular4tst.c: Modified Eular Method  連立オイラー法
 */
#include <stdio.h>
#include <math.h>
#define N 2

int Eular(int (*func)(double x[], double t, double f[]),
 double *t, double *x, double h, int n)
/*
    Differencial equation
    by using Eular method.
    Result is stored in "x".
 */
{
    int i;
    double f[N];

    (*func)(x, *t, f);
    for (i=0;i<n;i++) x[i] += h*f[i];
    *t += h;
    return 0;
}

double alpha, beta;
int rhs(double x[], double t, double f[])
/*  dx/dt = (f1(x,t),f2(x,t)) = (-a0*x0-a1*x1, -a0*x1-a1*x0) */
{
    double a0, a1;

    a0 = (alpha+beta)/2;
    a1 = (alpha-beta)/2;
    f[0] = -a0*x[0]-a1*x[1];
    f[1] = -a0*x[1]-a1*x[0];
    return 0;
}

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

//  Parameters
    alpha = 100; beta = 0.1;
    h = 0.01; tmax = 10;
    nprnt = 10;
    nstep = (int)(tmax/h+0.5); nskip = nstep/nprnt;
//  Define initial conditions
    t = 0; x[0] = 1; x[1] = 0;
//  Step on
    for (jstep=0; jstep<=nstep; jstep++)
    {
        if (jstep == (jstep/nskip)*nskip)
        {
            if (t < 700) a = exp(-alpha*t); else a = 0;
            b = exp(-beta*t);
            error[0] = x[0]-(a+b)/2;
            error[1] = x[1]-(a-b)/2;
            printf("%6.3f",t);
            for (i=0;i<N;i++) printf(" %14.7e %14.7e", x[i], error[i]);
            printf("\n");
        }
    //  --------------------------
        Eular(rhs, &t, x, h, N);
    //  --------------------------
    }
    return 0;
}

/*  ---------   実行例: 教科書 p.162 例4   ---------  */
/*
Z:\src\C>Eular4tst
 0.000  1.0000000e+00  0.0000000e+00  0.0000000e+00  0.0000000e+00
 1.000  4.5239607e-01 -2.2635461e-05 -4.5239607e-01  2.2635461e-05
 2.000  4.0932441e-01 -4.0961800e-05 -4.0932441e-01  4.0961800e-05
 3.000  3.7035352e-01 -5.5594263e-05 -3.7035352e-01  5.5594263e-05
 4.000  3.3509295e-01 -6.7070014e-05 -3.3509295e-01  6.7070014e-05
 5.000  3.0318947e-01 -7.5857426e-05 -3.0318947e-01  7.5857426e-05
 6.000  2.7432345e-01 -8.2364304e-05 -2.7432345e-01  8.2364304e-05
 7.000  2.4820571e-01 -8.6945180e-05 -2.4820571e-01  8.6945180e-05
 8.000  2.2457457e-01 -8.9907754e-05 -2.2457457e-01  8.9907754e-05
 9.000  2.0319331e-01 -9.1518598e-05 -2.0319331e-01  9.1518598e-05
10.000  1.8384771e-01 -9.2008200e-05 -1.8384771e-01  9.2008200e-05
*/