#include <iostream> 
#include <cmath>
using namespace std;

double A; // size where V acts
double V0; // height potential
double Ecur; // current trial value for the eigen energy of the wave function

double V(double x) { // infinite square well potential with internal barrier
	if (fabs(x)> 1.000001) return 1e50;
	if (fabs(x)<A) return V0;
	return 0;
}

/*
double V(double x) {
	return x*x*50; // H.O. Potential (for r=0.1x or for c=100). Eigenstates 10(n+1/2)
	return 0;
}
*/

double f(double x) { // returns nabla-squared = -2T = 2(V-E)
  return 2*(V(x)-Ecur);
}

int main() {
	double DELTA_X = 2/200.0; // choose stepsize such that it ends at the boundary of the square well
	A = 0.405; // choose size for barrier
	V0=10; // choose potential for barrier
	double Estart=1; // start value, first trial
	Ecur=Estart;
	double DeltaE = 1; // start step size in energies.
	double prevend; // value at x=1 of previous trial
	double psi0,psi1,psi2; // psi at k-1,k,k+1
	while (fabs(DeltaE)>1e-10) { // continue until energy is correct within 1e-10
		psi0=0;
		psi1=1; // arbitrary first value for psi. note, that the schroedinger equation is linear in psi. if psi is a solution, so is 100psi or -2psi. we calculate an unnormalized wave function
		double x =-1+ DELTA_X; // x at psi_k
		cout << Ecur << " " << DeltaE << " -1 0" << endl;
		cout << Ecur << " " << DeltaE << " " << x << " " << psi1 << endl;
		while (x<1.0-DELTA_X/2.0) {
			double fac = DELTA_X*DELTA_X/12;
			psi2=psi1*(2 + 10.0* fac*f(x)) -psi0*(1-fac*f(x-DELTA_X));
			x+=DELTA_X; // x at psi_k+1 = psi2
			psi2=psi2/(1-fac*f(x));
			cout << Ecur << " " << DeltaE << " " << x << " " << psi2 << endl;
			psi0=psi1;
			psi1=psi2;
		}
		if (Ecur==Estart) {
			prevend=psi2;
			Ecur+=DeltaE;
		} else {
			double df = psi2/(prevend-psi2)*DeltaE; // Newton-Rhapson 
			Ecur +=df;
			DeltaE = df;
			prevend=psi2;
		}
	}
	cout << endl;
	return 0;
}