T = 200;
h = 0.0005;
N = 40000;
QtMatrix = zeros(2,N+1);
PtMatrix = zeros(2,N+1);
%N = 40000;
E = 0.6;
Qprime = @(t) 1*PtMatrix(:,t);
%
%Pprime = @(tstep) (-(1/((QtMatrix(1,tstep)^2+ QtMatrix(2,tstep)^2)^(3/2))))*QtMatrix(:,tstep);
%Pprime = @(tstep) (-(1/((QtMatrix(tstep,1)^2+ QtMatrix(tstep,2)^2)^(3/2))))*QtMatrix(:,tstep);
Pprime = @(s) -(1/((QtMatrix(1,s)^2+QtMatrix(2,s)^2)^(3/2)))*QtMatrix(:,s);
Angular = @(ts) QtMatrix(ts,1)*PtMatrix(ts,2) + QtMatrix(ts,2)*PtMatrix(ts,1);
Hamilton = @(tst) ((PtMatrix(tst,1)^2+PtMatrix(tst,2)^2)/2)-(1/(sqrt((QtMatrix(tst,1)^2)+(QtMatrix(tst,2))^2)));
QtMatrix(1,1) = 1-E;
QtMatrix(2,1) = 0;
PtMatrix(1,1) = 0;
PtMatrix(2,1) = sqrt(1+E/1-E);
t = 0;
Avector = zeros(N+1,1);
Hvector = zeros(N+1,1);

%Euler's method part
for i=1:N
   % t = 0+i*h;
   % Ql = Qprime(i);
   %QtMatrix(:,i+1) = QtMatrix(:,i) + h*Qprime(i);
   QtMatrix(:,i+1) = QtMatrix(:,i) + h*1*PtMatrix(:,i);
   PtMatrix(:,i+1) = PtMatrix(:,i) + h*-(1/((QtMatrix(1,i)^2+QtMatrix(2,i)^2)^(3/2)))*QtMatrix(:,i);
end

figure;
scatter(PtMatrix(1,:), PtMatrix(2,:),5, 'r', '.')
hold('on');
scatter(QtMatrix(1,:),QtMatrix(2,:),5,'b','.')

% for z=1:N+1
%     Avector(z) = Angular(z);
%     Hvector(z) = Hamilton(z);
% end

for z=1:N+1
    Avector(z) = QtMatrix(1,z)*PtMatrix(2,z) + QtMatrix(2,z)*PtMatrix(1,z);
    Hvector(z) = ((PtMatrix(1,z)^2+PtMatrix(2,z)^2)/2)-(1/(sqrt((QtMatrix(1,z)^2)+(QtMatrix(2,z))^2)));
end