function [t,y] = pr2
%
% Solves the Perturbed Kepler problem with high accuracy. 
% Returns the solution on an equidistant grid. 
%
y0 = [0.4,0,0,2]
h = 0.1; 
tint = [0:h:500];

options = odeset('AbsTol',1.e-10,'RelTol',1.e-10);
[t,y] = ode45(@f,tint,y0,options);

plot(y(:,1),y(:,2));
grid
axis('square');
xlabel('q_1');
ylabel('q_2');

function yp = f(t,y)
% The perturbed Kepler problem

r = sqrt(y(1)^2+y(2)^2);
alpha = 0.01;
yp = zeros(4,1);
yp(1) = y(3);
yp(2) = y(4);
yp(3) = -y(1)/r^3-1.5*alpha*y(1)/r^5;
yp(4) = -y(2)/r^3-1.5*alpha*y(2)/r^5;

    