% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -
% Vers. 4.0 2016
%%
% ausgewaehlte Startwerte
y0(2)=0;
% Zufallswerte fuer die Impulse
zwi=2*rand;
z(randperm(2))=[zwi;sqrt(4-zwi^2)*rand];
y0(3)=z(1);
y0(4)=z(2);
epsilon=-0.5;   % sinnvolle Werte sind negativ.
% Aus der Energieerhaltung folgt
y0(1)=sqrt((1/2*(y0(3)^2+y0(4)^2)-2)/epsilon);
options = odeset('RelTol',1e-11,'AbsTol',1e-11,'MaxStep',0.1);
tmax=200;
[t,y] = ode23(@(t,y) diahdgl(t,y,epsilon),[0 tmax],y0,options);