%% ENM 7.13

U0 = [1 0 0 0.5];
T = 3;

dt = 0.000001;
n_steps = T/dt;
U = zeros(4,n_steps+1);
U(:,1) = U0;
for i = 1:n_steps
    r = sqrt(U(1,i)^2 + U(3,i)^2);
 U(:,i+1) = U(:,i) + dt * [U(2,i); -U(1,i)/(r^3); U(4,i); -U(3,i)/(r^3)];

end
figure
plot(0,0,'*')
hold on
axis equal
plot(U(1,:),U(3,:))
