clear all                           % Programm zur numerischen Berechnung der zeitabhängigen Wellenfunktion eines Atoms mit Z Elektronen
clc                                 % Programm funktioniert mit beliebiger Anzahl an Elektronen

Z = 1;                              % Ordnungs-/Elektronen-/Protonenzahl des Atomes
e0 = 8.85418781762*10^(0-12);       % Permittivität
me = 9.10938291*10^(0-31);          % Elektronenmasse
h = 6.62606957*10^(0-34);           % planksches Wirkungsquantum
e = 1.602176565*10^(0-19);          % Elementarladung
a0 = 0.52917721092*10^(0-10);       % bohrscher Radius
I = [-1:0.02:1].*(a0.*4);           % Intervall zur numerischen Auswertung Differential
I0 = I(2:(end-1));                  % Intervall zur numerischen Auswertung 2 Elemente weniger, da 2 beim Differenzenquozienten verloren gehen x
I_1= I(1:(end-2));                  % Intervall zur numerischen Auswertung 2 Elemente weniger, da 2 beim Differenzenquozienten verloren gehen x-dx
I1 = I(3:end);                      % Intervall zur numerischen Auswertung 2 Elemente weniger, da 2 beim Differenzenquozienten verloren gehen x+dx
Ik = [I1;I0;I_1];                   % alle 3 Intervalle werden in einer Matrix abgespeichert 
t = [0:0.01:1];                     % Intervall Zeit
dIndizes = 2:(numel(I)-1);          % Vektor zur Indizierung für nicht betroffene x,y,z 1.Ableitung
ddIndizes = 1:(numel(I)-2);         % Vektor zur Indizierung für nicht betroffene x,y,z 2.Ableitung

for k=1:(3.*Z)
   
    dimension = ones(1,3.*Z);              % für reshape, richtet das Intervall in die gewünschte Richtung (Dimension) aus
    dimension(k) = numel(I0);
    matrix = ones(1,3.*Z).*numel(I0);       % für repmat, macht aus dem Vektor eine Matrix in 3Z Dimensionen
    matrix(k) = 1;
    
    X0 = repmat(reshape(I0,dimension),matrix); % 3Z dimensionale Matrix für jede Richtung; enthält Intervall für phi("x")
   
    matrizen0{k} = X0;                        % abspeichern der Intervallmatrizen in eine Zelle, zum leichteren auslesen
  
end

for k=1:(3.*Z)
    
    matrizen=matrizen0;
    
    for m=1:3
    
        dimension = ones(1,3.*Z);               % für reshape, richtet das Intervall in die gewünschte Richtung (Dimension) aus
        dimension(k) = numel(I0);
        matrix = ones(1,3.*Z).*numel(I0);       % für repmat, macht aus dem Vektor eine Matrix in 3Z Dimensionen
        matrix(k) = 1;
    
        X0 = repmat(reshape(Ik(m,:),dimension),matrix); % 3Z dimensionale Matrix für jede Richtung; enthält Intervall für phi("x+dx/x/x-dx"), Auswahl erfolgt mit m 
    
        matrizen{k}=X0;
    
        phi = exp(0-sqrt(matrizen{1}.^2+matrizen{2}.^2+matrizen{3}.^2.)/a0); % Bsp. für H; allgemein : xk=3*(k-1)+1, yk=3*(k-1)+2, zk=3*k

        if m==2                     % Bestimmung des Vorfaktors
            a=0-2;
        else
            a=1;
        end
        
        if k==1 && m==2             % Abspeicherung von phi0
            phi0=phi;
        end
        
        if k==1 && m==1
           ddphi_sum=zeros(size(phi)); 
        end
        
        ddphi_sum=ddphi_sum+(a.*phi)./(0.02.*4.*a0).^2;     % Summenbildung über d²phi/dx² = (phi(x+dx)-2*phi(x)+phi(x-dx))/(Schrittweite des Intervalles)²
        
    end
end

for k=1:Z                                % Berechnung der Abstände Elektron-Kern + Summenbildung aller Terme
    
    rNE = (matrizen0{((k-1).*3+1)}.^2+matrizen0{((k-1).*3+2)}.^2+matrizen0{((k).*3)}.^2).^(1./2);
    
    if k==1
        rNE_sum = zeros(size(rNE));
    end
    
    rNE_sum = rNE_sum + 1./rNE;
    
    if Z==1                                  % Für Z=1 existiert keine Elektron-Elektron- Wechselwirkung --> 0
    rEE_sum = zeros(size(rNE_sum));
    end
    
end

for k=1:(Z-1)                            % Berechnung der Abstände Elektron-Elektron + Summenbildung aller Terme
    for m=2:Z
        
        rEE = ((matrizen0{((k-1).*3+1)}-matrizen0{((m-1).*3+1)}).^2+(matrizen0{((k-1).*3+2)}-matrizen0{((m-1).*3+2)}).^2+(matrizen0{((k).*3)}-matrizen0{((m).*3)}).^2).^(1./2);
        
        if k==1 && m==2
            rEE_sum = zeros(size(rEE));
        end
        
        rEE_sum = rEE_sum + 1./rEE;  
        
    end
end

phi_dt = (rEE_sum.*(e.^2./(h.*2.*e0)).*phi0)./(1i) - (rNE_sum.*(Z.*e.^2./(h.*2.*e0)).*phi0)./(1i) - ((h./(4.*pi.*me)).*ddphi_sum)./(1i); % Wellenfunktion d(phi)/dt

for k=1:numel(t)
    
    PHI_t{k} = phi_dt.*t(k);             % Lösung der zeitabhängigen Wellenfunktion; die Wellenfunktion zum Zeitpunkt t(k) ist in PHI_t{k} gespeichert
    
end