#include "stdafx.h"
#include <iostream>
#include <stdio.h>
#include <cstdlib>
#include <cmath>

int N = 3;
int count = 20;
double y(double x) {return N*exp(x);}
double f(double x) {return 21*exp(x);}
double p(double x) {return N;}
double q(double x) {return N;}
double* TridiagonalMatrixAlgorithm(double* a, double* b, double* c, double* d) {
	int n = count-2;
	double *alpha = new double[n-1];
	double *beta = new double[n-1];
	double *x = new double[n];
	alpha[0] = (-1.0) * c[0] / b[0];
	beta[0] = d[0] / b[0];

	for (int i = 1; i < n - 1; i++) 	{
		alpha[i] = ((-1.0) * c[i]) / (a[i - 1] * alpha[i - 1] + b[i]);
		beta[i] = (d[i] - a[i - 1] * beta[i - 1]) / (a[i - 1] * alpha[i - 1] + b[i]);
	}
	x[n - 1] = (d[n - 1] - beta[n - 2] * a[n - 2]) / (a[n - 2] * alpha[n - 2] + b[n - 1]);
	for (int i = n - 2; i >= 0; i--) {
		x[i] = alpha[i] * x[i + 1] + beta[i];
	}
	return x;
}

int main(int argc, char* argv[]) {
	
	double *res = new double[count];
	double *res1 = new double[count];
	double *x = new double[count];
	double l = 0;
	double r = 1;
	double h = (r-l)/(double)(count-1);
	for (int i = 0; i < count; i++) {
		x[i] = h * i;
    }

	res[0] = y(x[0]);
	res[count-1] = y(x[count-1]);

	double *a = new double[count-3];
	double *b = new double[count-2];
	double *c = new double[count-3];
	double *d = new double[count-2];

	b[0] = q(x[1])*h*h-2.0;
	c[0] = (1.0+p(x[1])*h/2.0);
	d[0] = f(x[1])*h*h - (1-p(x[1])*h/2.0)*res[0];
	for (int i=1;i<count-3;i++) {
		a[i-1] = 1-p(x[i+1])*h/2.0;
		c[i] = 1+p(x[i+1])*h/2.0;
		b[i] = q(x[i+1])*h*h-2.0;
		d[i] = f(x[i+1])*h*h;
	}

	a[count-4] = 1-p(x[count-2])*h/2.0;
	b[count-3] = q(x[count-2])*h*h-2.0;
	d[count-3] = f(x[count-2])*h*h - (1-p(x[count-2])*h/2.0)*res[count-1];

	res1 = TridiagonalMatrixAlgorithm(a, b, c, d);
	for (int i=0;i<count-2;i++) {
		res[i+1] = res1[i];
	}

	for (int i=0;i<count;i++) {
		printf("%f\t, %f\t, %f\t, %f\n", x[i], y(x[i]), res[i], abs(y(x[i])-res[i]));
	}

}