function twobody

% This is a Matlab simulation of the motion in a
% concervative central force field
%
% X''= - F(X)= -F(|X|)X = - X/|X|^p
% 
% Copyleft 2009 - Radu C. Cascaval, UCCS

clf
X0=[1;0.99;1;0];
opts=odeset('reltol', 1.e-6); % sets the relative tolerance for ode45
[t,X]=ode45(@(t,X) odefun(X), [0,100], X0,opts); 

plot(0,0,'or','MarkerSize',11), hold on
quiver(X0(1),X0(3),X0(2),X0(4),'Color',[0,0.5,0],'LineWidth',5');
F=odefun(X0);
quiver(X0(1),X0(3),F(2),F(4),'Color',[0.5,0,0],'LineWidth',5');

plot(X(:,1),X(:,3),'-');
axis equal
xlabel('x1')
ylabel('x2')
grid on
hold off
end

function dXdt=odefun(X)

% You may try any of the following values of p

p = 3; % classical Newtonian force 
% p = 3.2; % non-Newtonian central force
% p = 4; % non-Newtonian central force - black hole

r=sqrt(X(1)^2+X(3)^2);

dXdt=[X(2); -X(1)/r^p; X(4); -X(3)/r^p];
end