% Grundenergie der Elektronen des Helium Atoms
% Golubkov Andrej

n=1600;
rmax=6;
h=rmax/n;

i=[1:n];
ii=[1:n-1];
r=i*h;

u=sqrt(32) * r .* exp(-2*r).*2;

%spline für u^2
pp = csape(r,u.^2,'variational');
spline_coeff = pp.coefs;
a1=spline_coeff(:,4);
b1=spline_coeff(:,3);
c1=spline_coeff(:,2);
d1=spline_coeff(:,1);

% u^2 normierung
I = sum(a1*h + 0.5*b1*h^2 + 1/3*c1*h^3 + 1/4*d1*h^4)
u = u/sqrt(I);

%neue spline für u^2
pp = csape(r,u.^2,'variational');
spline_coeff = pp.coefs;
a1=spline_coeff(:,4);
b1=spline_coeff(:,3);
c1=spline_coeff(:,2);
d1=spline_coeff(:,1);

%spline für u^2/r
pp = csape(r,u.^2./r,'variational');
spline_coeff = pp.coefs;
A1=spline_coeff(:,4);
B1=spline_coeff(:,3);
C1=spline_coeff(:,2);
D1=spline_coeff(:,1);


% I2(g) Berechnen mit g = 1...n, I2(1) = 0, I2(n)=int(u^2,0,inf)
I2_einzel(1) = 0;
I2_einzel(ii+1) = a1(ii)*h + 0.5*b1(ii)*h^2 + 1/3*c1(ii)*h^3 + 1/4*d1(ii)*h^4;
I2 = cumsum(I2_einzel);



% I3(g) Berechnen mit g = 1...n, I3(n) = 0, I3(1)=int(u^2/r,0,inf)
I3_einzel(1) = 0;
I3_einzel(ii+1) = A1(ii)*h + 0.5*B1(ii)*h^2 + 1/3*C1(ii)*h^3 + 1/4*D1(ii)*h^4;
I3_einzel = fliplr(I3_einzel); % Spiegeln
I3 = cumsum(I3_einzel);
I3 = fliplr(I3);




