#include <iostream>
#include <iomanip>
#include <math.h>

using namespace std;

class STMatrix {
private:
	enum {normal, upper, lower} state;
	int n;
	double* stm;
public:
	STMatrix(int n) {
		this->n = n;
		stm = new double[n*2-1];
		for (int i=0; i<n*2-1; i++)
			stm[i] = 0;
		state = normal;
	}
	~STMatrix() {
		delete stm;	
	}
	int getN() {
		return n;
	}
	double getElement(int row, int col) {
		if (state == lower && col>row)
			return 0;
		if (state == upper && col<row)
			return 0;
		if (col == row+1) {
			col--; //adjust for symmetricity
			row++;
		}
		if (col < row-1 || col > row+1) {
			return 0;
		}
		if (col == row)
			return stm[row];
		else
			return stm[n+col];
	}
	void setElement(int row, int col, double ele) {
		if (col == row+1) {
			col--; //adjust for symmetricity
			row++;
		}
		if (col < row-1 || col > row+1) {
			cerr << "trying to set an off-diagonal element" << endl;
			return; //can't happen, but if it does let's not crash
		}
		if (col == row)
			stm[row] = ele;
		else
			stm[n+col] = ele;
	}
	void eliminate() {
		for (int k=0; k<n-1; k++) {
			double kpkp = getElement(k+1, k+1);
			double kkp = getElement(k, k+1);
			double kk = getElement(k, k);
			setElement(k+1, k+1, kpkp-kkp*kkp/kk);
		}
		state = upper;
	}
};

double potential(double x, int which) {
	if (which == 0)
		return x*x;
	return 2*pow(M_E, x);	
}

double solution(double x, int which) {
	if (which == 0)
		return x*x+2;
	return (1-x)*pow(M_E, x);	
}

void productST(STMatrix& m, double* x, double* y) {
	for (int i=0; i<m.getN(); i++) {
		y[i] = 0;
		for (int j=i-1; j<i+2; j++) {
			if (j>=0 && j<m.getN()) {
				y[i] += x[j]*m.getElement(i, j);
			}
		}
	}
} 

void solveST(STMatrix& m, double* b) {
	m.eliminate();
	int n = m.getN();
	b[0] /= m.getElement(0, 0);
	for (int i=1; i<n; i++)
		b[i] = (b[i]-m.getElement(i-1, i)*b[i-1])/m.getElement(i, i);
	for (int i=0; i<n; i++)
		b[i] *= m.getElement(i, i);
	b[n-1] /= m.getElement(n-1, n-1);
	for (int i=n-2; i>=0; i--)
		b[i] = (b[i]-m.getElement(i, i+1)*b[i+1])/m.getElement(i, i);
}

double BVP(int k, int which, double alpha, double beta) {
	int n = (int)pow(2, k);
	double h = 1/((double)n+1);
	STMatrix m(n);
	double* b = new double[n];
	for (int i=0; i<n; i++) {
		m.setElement(i, i, 2+h*h);
		if (i == n-1)
			m.setElement(i, i-1, -1);
		else
			m.setElement(i, i+1, -1);
		b[i] = h*h*potential((i+1)*h, which);
	}
	b[0] += alpha;
	b[n-1] += beta;
	
	solveST(m, b);
	
	double maxdiff = 0;
	for (int i=0; i<n; i++) {
		double diff = solution(i*h, which)-b[i];
		if (diff<0)
			diff *= -1;
		if (diff>maxdiff)
			maxdiff = diff;
	}
	
	delete b;
	
	return maxdiff;
}

int main(int argc, char** argv) {
	cout << setprecision(15);
	cout << "a)\tk\t||e||_inf" << endl;
	for (int k=1; k<=16; k++)
		cout << "\t" << k << "\t" << BVP(k, 0, 2, 3) << endl;
	cout << endl;
		
	cout << "b)\tk\t||e||_inf" << endl;
	for (int k=1; k<=16; k++)
		cout << "\t" << k << "\t" << BVP(k, 1, 1, 0) << endl;
	
	return 0;	
}