%% Struktur des Programms Elektrik.m
Tend=30;h=0.12;N=floor(Tend/h);
% Anfangswerten q0,p0 + Anfangsenergie ham0
q0=zeros(3,1);q0(1)=0.2;q0(2)=0.2;q0(3)=sqrt(0.92);
p0=zeros(3,1);p0(1)=1;p0(2)=-1;p0(3)=0;
ham0=...
% Erste Approximation mit explitizes Euler
qit=zeros(3,1);pit=zeros(3,1);Qit=zeros(3,1);Pit=zeros(3,1);
qit=q0+h.*[p0(1)+q0(2);p0(2)-q0(1);p0(3)];
pit=p0;
Qit=q0+.5*h.*[p0(1)+q0(2);p0(2)-q0(1);p0(3)];
Pit=p0;
lambda0it=1.439;lambda1it=lambda0it;

F=zeros(14,1);Fprime=zeros(14,14);
y0=[Qit;Pit;qit;lambda0it;pit;lambda1it];

for j=1:N
  err1=1;err2=1;
  while ((abs(err1) > 10^(-14)) & (abs(err2) > 10^(-14)))
    %%%%%%%%%%%%%%%%%%%%
    %% (1,1)-GL SPARK %%   
    %%%%%%%%%%%%%%%%%%%%
    F(1:3)=.5.*(qit+q0)-Qit;
    F(4)=p0(1)-.5*h*(Qit(1)-Pit(2))- ... 
         h*q0(1)*lambda0it-Pit(1);
    F(5)=p0(2)-.5*h*(Qit(2)+Pit(1))- ... 
         h*q0(2)*lambda0it-Pit(2);
    F(6)=p0(3)+.5*h-h*q0(3)*lambda0it-Pit(3);
    F(7)=q0(1)+h*(Pit(1)+Qit(2))-qit(1);
    F(8)=q0(2)+h*(Pit(2)-Qit(1))-qit(2);
    F(9)=q0(3)+h*Pit(3)-qit(3);
    F(10)=qit(1)^2+qit(2)^2+qit(3)^2-1;
    F(11)=p0(1)-h*(Qit(1)-Pit(2))-h*q0(1)*lambda0it ... 
         -h*qit(1)*lambda1it-pit(1);
    F(12)=p0(2)-h*(Qit(2)+Pit(1))-h*q0(2)*lambda0it ... 
         -h*qit(2)*lambda1it-pit(2);
    F(13)=p0(3)+h-h*q0(3)*lambda0it ... 
         -h*qit(3)*lambda1it-pit(3);
    F(14)=2*(qit(1)*(pit(1)+qit(2))+ ... 
          qit(2)*(pit(2)-qit(1))+ ... 
          qit(3)*pit(3));
        
    Fprime(1,:)=[-1 0 0 0 0 0 1/2 0 0 0 0 0 0 0];
    Fprime(2,:)=[0 -1 0 0 0 0 0 1/2 0 0 0 0 0 0];
    Fprime(3,:)=[0 0 -1 0 0 0 0 0 1/2 0 0 0 0 0];
    Fprime(4,:)=[-h/2 0 0 -1 h/2 0 0 0 0 -h*q0(1) ... 
                 0 0 0 0];
    Fprime(5,:)=[0 -h/2 0 -h/2 -1 0 0 0 0 -h*q0(2) ... 
                 0 0 0 0];
    Fprime(6,:)=[0 0 0 0 0 -1 0 0 0 -h*q0(3) ... 
                 0 0 0 0];
    Fprime(7,:)=[0 h 0 h 0 0 -1 0 0 0 0 0 0 0];
    Fprime(8,:)=[-h 0 0 0 h 0 0 -1 0 0 0 0 0 0];
    Fprime(9,:)=[0 0 0 0 0 h 0 0 -1 0 0 0 0 0];
    Fprime(10,:)=[0 0 0 0 0 0 2*qit(1) 2*qit(2) ... 
                  2*qit(3) 0 0 0 0 0];
    Fprime(11,:)=[-h 0 0 0 h 0 -h*lambda1it 0 0 ... 
                  -h*q0(1) -1 0 0 -h*qit(1)];
    Fprime(12,:)=[0 -h 0 -h 0 0 0 -h*lambda1it 0 ... 
                  -h*q0(2) 0 -1 0 -h*qit(2)];
    Fprime(13,:)=[0 0 0 0 0 0 0 0 -h*lambda1it ... 
                  -h*q0(3) 0 0 -1 -h*qit(3)];
    Fprime(14,:)=[0 0 0 0 0 0 ... 
                  2*(pit(1)+qit(2))-2*qit(2) ... 
		  2*qit(1)+2*(pit(2)-qit(1)) ... 
		  2*pit(3) 0 ... 
		  2*qit(1) 2*qit(2) 2*qit(3) 0];
    Fprime=sparse(Fprime);
    % Newton Verfahren + Berechnung des Fehlers
    ...	  
  end  
  % Berechnung der Energie+Einspeichern für die Plots
  % Update die Daten
  ...
end
% Plot der Energie+Visulaisierung
...;fig=figure;set(fig,'doublebuffer','on')
xlim([-1 1]);ylim([-1 1]);zlim([-1 1])

p=plot3(q1plot(1),q2plot(1),q3plot(1),'ro');
for j=2:N
  set(p,'xdata',q1plot(1:j),'ydata',q2plot(1:j),'zdata',q3plot(1:j));
  drawnow;
end  
