%% Bahn parallel zum Magnetfeld

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -
% Vers. 1.0 2016

% ausgewaehlte Startwerte
y0=[0.0001,0.];
% 
epsilona=linspace(-0.5,-0.09,200);   % sinnvolle Werte sind negativ.
% Aus der Energieerhaltung folgt
% options = odeset('RelTol',1e-11,'AbsTol',1e-11,'MaxStep',0.1);
for n=1:length(epsilona)
    n
    epsilon = epsilona(n);
    tmax=4*pi/sqrt(-2*epsilon);
    options = odeset('RelTol',1e-10,'AbsTol',1e-10,'MaxStep',0.1,'Events',@(t,y) diahparallelevent(t,y,epsilon));
    [t,y] = ode45(@(t,y) diahparalleldgl(t,y,epsilon),[0 tmax],y0,options);
    ya = y(:,1);
    ya = ya(t>t(end)/2);
    yv(n)=max(ya)/y(1,1);
    clear t y
end
plot(epsilona,yv),shg
yv(yv>1) = 2;
yv(yv<=1) = 0;
figure
stairs(epsilona,yv);ylim([0,2.05]),shg