Code

#include <fcntl.h>
#include <math.h>
#include <omp.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>

#define F(y) (exp(-(y)))
#define dF(y) (-exp(-(y)))
#define X(i) ((x0) + (i)*(h))

void reduce(int N, double h, double *y, double *A, double *B, double *C,
            double *D, double *d)
{
    double p = h * h / 12.0;
    double q = h * h * 5.0 / 6.0;
    int stride = 1;

    A[0] = 0; B[0] = 1.0; C[0] = 0; D[0] = 0;

    A[N-1] = 0; B[N-1] = 1.0; C[N-1] = 0; D[N-1] = 0;

    #pragma omp parallel for
    for (int i = 1; i < N - 1; i++) {
        A[i] = 1 - dF(y[i - 1]) * p;
        B[i] = -2 - dF(y[i]) * q;
        C[i] = 1 - dF(y[i + 1]) * p;
        D[i] = (y[i - 1] * dF(y[i - 1]) + F(y[i - 1])) * p +
               (y[i] * dF(y[i]) + F(y[i])) * q +
               (y[i + 1] * dF(y[i + 1]) + F(y[i + 1])) * p -
               (y[i - 1] - 2 * y[i] + y[i + 1]);
    }

    for (int n = N, low = 2; n > 1; n >>= 2, low <<= 1, stride <<= 1) {
        #pragma omp parallel for
        for (int i = low - 1; i < N; i += stride * 2) {
            double alpha = -A[i] / B[i - stride];
            double beta  = -C[i] / B[i + stride];

            A[i] = alpha * A[i - stride];
            B[i] = B[i] + alpha * C[i - stride] + beta * A[i + stride];
            C[i] = beta * C[i + stride];
            D[i] = D[i] + alpha * D[i - stride] + beta * D[i + stride];
        }
    }

    d[0] = 0; d[N/2] = D[N / 2] / B[N / 2]; d[N - 1] = 0;

    for (stride >>= 1; stride >= 1; stride >>= 1) {
        #pragma omp parallel for
        for (int i = stride - 1; i < N; i += stride * 2)
            d[i] = (D[i] - ((i - stride > 0) ? (A[i] * d[i - stride]) : 0.)
                    - ((i + stride < N) ? (C[i] * d[i + stride]) : 0.)) / B[i];
    }
}

int main(int argc, char** argv)
{
    double b0 = 0, b1 = 1.7, db = 0.1;
    double x0 = 0, xN = 1.0;
    double eps = 0.000001;
    int N = 3000;

    double h = fabs(xN - x0) / N;

    double *A = malloc((N + 1) * sizeof(double));
    double *B = malloc((N + 1) * sizeof(double));
    double *C = malloc((N + 1) * sizeof(double));
    double *D = malloc((N + 1) * sizeof(double));
    double *d = malloc((N + 1) * sizeof(double));
    double *y = malloc((N + 1) * sizeof(double));

    int out_file = open("Kaziakhmedov.txt", O_WRONLY | O_CREAT);

    dprintf(out_file, "plot ");
    /*
     * Draw parameter range
     */
    for (double b = b0; b < b1; b += db)
        dprintf(out_file, "'-' pt 7 ps 0.25 title 'b = %2.2lf', ", b);
    dprintf(out_file, "'-' pt 7 ps 0.25 title 'b = %2.2lf'\n", b1);
    /*
     * Main calc
     */
    for (double b = b0; b <= b1 + eps; b += db) {
        #pragma omp parallel for
        for (int i = 0; i < N + 1; i++)
            y[i] = 1.0 + X(i) * (b - 1.0);

        while (1) {
            double max_delta = -1.;

            reduce(N + 1, h, y, A, B, C, D, d);
            #pragma omp parallel for reduction(max:max_delta)
            for (int i = 1; i < N; i++) {
                y[i] += d[i];
                if (fabs(d[i]) > max_delta)
                    max_delta = fabs(d[i]);
            }
            if (max_delta < eps)
                break;
        }

        for (int i = 0; i < N + 1; i++)
            dprintf(out_file, "%lf %lf\n", X(i), y[i]);

        dprintf(out_file, "e\n");
    }

    close(out_file);
    return 0;
};