function dy = keplerDGLART(t,y,m1,Vorf,m2,G)
% Differentialgleichung

% Vers.2.0  21.12.2013  --> ART-Korrektur 
% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -


[y2art,y3art] = ARTkorr(y,m1,m2,G);


dy(1) = y(3)/m1;
dy(2) = y(4)/(m1*y(1)^2) + y2art;
dy(3) = y(4)^2/(m1*y(1)^3) - Vorf/(y(1)^2) + y3art;
dy(4) = 0;
dy = dy.';

function [y2art,y3art] = ARTkorr(y,m1,m2,G)

% N??herungsweise Berechnung der Korrekturen der allg. Relativit??tstheorie

cq = 2.99792e08^2; % Lichgeschwindigkeit in m/s zum Quadrat
mf = G*((m1+m2)^2)/(m1*m2*cq*y(1)^3); 

y2art = -2*y(4)*mf;
y3art = -3*y(4)^2*mf/y(1);
