% This is the code for Problem 7 of HW 2.
% Comments in  lower-case  letters are provided for your information only.
% You need to follow the directions typed in CAPITAL LETTERS.

clear all
t0=0;                         % initial time
tf=100;                        % final time
                              %%%--->>> YOU MAY NEED TO VARY THIS NUMBER TO OBTAIN 
                              %%%--->>>  A CLOSED ORBIT (OR A PIECE OF THE ORBIT
                              %%%--->>>  CLEARLY SHOWING THAT IT IS NOT CLOSED)
initial_position=[0, 1];      % These are the initial x- and y-coordinates of the planet.
                              %  The Sun is assumed to be located at (0,0).
vy0=input('vy0 = ')
initial_velocity=[-1, vy0];    % These are the initial x- and y-velocities of the planet.
                              %%%--->>> YOU NEED TO SUPPLY A NUMERIC VALUE FOR vy0.
                              %
                              %  Note that the initial x-velocity should not equal 0,
                              %  for if it is, then the planet would move either 
                              %  straight towards or straight away from the Sun, since
                              %  x(0)=0.
                              %  Then, as long as  vx0 is nonzero, one can set it to 1
                              %  without loss of generality. Indeed, setting it to another
                              %  value will simply rescale the time in the equation of motion.
                              %  (The previous sentence assumes that the initial condition 
                              %   is fixed. Changing _both_ the initial y-coordinate and
                              %   x-velocity will in general change the orbit, not just
                              %   rescale the time.)
z0=[initial_position, initial_velocity];
                              % this is a bookkeeping step
options=odeset('RelTol',10^(-9));
                              % increase the relative accuracy of ode45 (see below)
[t,z]=ode45('hw2_p7_funorbit',[t0 tf], z0,options);
                              % this solves the equation indicated in hw2_p6_funorbit.m
                              %  using Matlab's built-in ODE solver  ode45.m 

figure(271);                              
plot(z(:,1),z(:,2))           % plot the orbit
xlabel('x','fontsize',18)
ylabel('y','fontsize',18)