/* prog7a.c -- Solve the two-point boundary value problem

       y''+P(x)y'+Q(x)y=R(x), y(x0)=y0, y(xn)=yn

   by the central difference method.
   Written May 1, 2001 by Eric Olson for Mathematics 483

   For details see Section 5.2 in Numerical Analysis for Applied
   Mathematics, Science and Engineering by Greenspan and Casulli. */

#include <stdio.h>
#include <math.h>

#define Nmax 1024

typedef double vector[Nmax];
typedef vector tridiag[3];

static double P(double x){
    return cos(x);
}

static double Q(double x){
    return -(x * x + 1);
}

static double R(double x){
    return 2.0;
}

static trisolve(vector *a, double *b, int n, double *y){
    int i;
    double p;

    for (i = 1; i <= n - 2; i++) {
        p = a[0][i] / a[1][i - 1];
        a[1][i] -= p * a[2][i - 1];
        b[i] -= p * b[i - 1];
    }
    for (i = n - 1; i >= 2; i--) {
        p = a[2][i - 2] / a[1][i - 1];
        b[i - 2] -= p * b[i - 1];
    }
    for (i = 0; i <= n - 2; i++)
        y[i] = b[i] / a[1][i];
}

static CD(double x0, double y0, double xn, double yn, int n, double *y){
    static char *fmt = "%16.10lg %20.10lg\n";
    tridiag a;
    vector b;
    double x, h, hh2;
    int i;

    h = (xn - x0) / n;
    hh2 = 2.0 * h * h;
    for (i = 0; i <= n - 2; i++) {
        x = x0 + (i + 1) * h;
        b[i] = hh2 * R(x);
        a[0][i] = 2 - h * P(x);
        a[1][i] = hh2 * Q(x) - 4.0;
        a[2][i] = 2 + h * P(x);
    }
    b[0] -= a[0][0] * y0;
    b[n - 2] -= a[2][n - 2] * yn;

    trisolve(a, b, n, y);

    printf("#xi yi\n");
    printf(fmt, x0, y0);
    for (i = 1; i < n; i++) {
        x = x0 + i * h;
        printf(fmt, x, y[i - 1]);
    }
    printf(fmt, xn, yn);
    putchar('\n');
}

main(){
    int i, n;
    vector y;
    double x0, y0, xn, yn;

    x0 = -2.0; y0 = -1.0;
    xn = 2.0; yn = -1.0;

    n = 4;
    for (i = 2; i <= 8; i++) {
        CD(x0, y0, xn, yn, n, y);
        n *= 2;
    }
}