clear;
close all;

startVal = [ ... %% ... stands for look at next line
1,3,4,7,8,10,11,12; ...
];

startVal = startVal - i*1e-6;

[numEigenValues,n] = size(startVal);
cmap = jet(numEigenValues);

figureManyParticles = 16;

g_end = 2; %Maximaler Wert für die Kopplung g
numSteps = 300; %Anzahl Iterationsschritte


for  p = 1 : numEigenValues %Falls mehrere Startwerte gleichzeitig ausgewertet werden

y = startVal(p,:);
    
z = NaN(numSteps,1);
z(1)=0;
Energy = sum(real(y));

%Schleife zur Lösung des NGLs. Wie in BetheBCS.m schon auch hier mit
%"erwzungener" komplexer Darstellung (insbes. Z. 33,35)
for k =1:numSteps
    g = (g_end/numSteps)*k;   
    z(k+1)=g;
    f =@(x)BetheBCS(x,g);
    v  = [real(y(k,:));imag(y(k,:))];
    [res,~]=fsolve(f,v);
    y(k+1,:) = res(1,:) + 1i*res(2,:);
    Energy(k+1) =sum(real(y(k+1,:)));
end
    
%Energie, Realteil und Imaginärteil plotten (v gegen g)
figure(figureManyParticles);
subplot(2,2,1)
plot(z,Energy,'.','Color',cmap(p,:));
xlabel('Parameter g');
ylabel('Energy [a.u.]');
subplot(2,2,2)
plot(z,real(y),'b');
xlabel('Parameter g');
ylabel('Real(v)');
subplot(2,2,3)
plot(z,imag(y),'b')
xlabel('Parameter g');
ylabel('Imag(v)');
   
% set(gca,'ytick',[]);

end