function dy = diahphidgl(t,y,epsilon)
% Differentialgleichung zu nichtverschwindendem Drehimpuls

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


dy(1) = y(5);
dy(2) = y(6);
dy(3) = y(7);
dy(4) = y(8);
%
z12 = y(1).^2 + y(2).^2;
z34 = y(3).^2 + y(4).^2;
z1234 = z12 + z34;
%
zwi12 = (2*epsilon - 1/4*z34*(z1234 + z12));
zwi34 = (2*epsilon - 1/4*z12*(z1234 + z34));
%
dy(5) = zwi12 * y(1);
dy(6) = zwi12 * y(2);
dy(7) = zwi34 * y(3);
dy(8) = zwi34 * y(4);

dy = dy';
%0 = y(2)*y(5) - y(1)*y(6) - y(4)*y(7) + y(3)*y(8);
