eulertest.m

N = 100:10:10000;
% creates a vector for N, the different values for the number of steps 

error = zeros(size(N));
% creates an array to store the error in
for i=1:length(N)
    S = planetsolve(1, 2, 1, N(i));  
    error(i) = abs(S(end) - 1);
end
% calculates modulus of R(2pi) - 1 for every value of N in vector above

plot(1./N, error, '*') % plots error as a function of 1/N
xlabel('1/N')
ylabel('|R(2*pi)-1|')   
title('Error in R as a function of 1/N')

% the error approximates a straight line as 1/N tends to 0
% this implies that the error in the Euler method is proportional to the
% step size, (so as we increase the number of steps, the error decreases; 
% the larger N is, the better the approximation)ul

planetsolve.m

function R = planetsolve(h, gam, V0, N) 
% creates a function that inputs h, gamma, V0 and N (number of steps) and
% outputs R

deltax = (2*pi)/N;
% deltax is the difference between each step/ size of each step
% 2*pi represents the maximum value that theta can be


rEuler = zeros(1, N+1); % empty vectors to store solutions for R and V
vEuler = zeros(1, N+1);
rEuler(1) = 1; vEuler(1) = V0; % satisfies initial conditions           

for i = 1:N
    rEuler(i+1) = rEuler(i) + deltax*(((rEuler(i)^2) *vEuler(i))/h);
    vEuler(i+1) = vEuler(i) + deltax*((h/rEuler(i))-(gam/h));
end
% for loop carries out Euler's method

R = rEuler; % outputs solution for R

solvescript.m

syms u(x) v(x) h gam V0  % x represents theta
% identifies variables as symbolic
S = dsolve(diff(u, x) == -v/h, diff(v, x) ==  h*u - gam/ h, u(0) ==1, v(0) == V0); 
A = simplify(S.u, 100);
B = simplify(S.v, 100);
disp(A)
disp(B)
% solves system of ODEs using initial conditions
% displays simplified solution

N = 100; % uses 100 steps
xrange= linspace(0, 2*pi, N+1); % sets up range of x for Euler solution
rEuler = planetsolve(1, 2, 1, N); % finds solution for R as a function of x

xExactRange = linspace(0, 2*pi, 1000); % range of x for exact solutionn
uExact = sym(zeros(1,1000)); % zero vector to store u in
a1 = subs(A, h, 1); a2 = subs(a1, gam, 2); a3 = subs(a2, V0, 1);
%substitues values for h, gamma and V0 into solution

for i = 1:1000
    uExact(i) = subs(a3, x, xExactRange(i));
end
rExact = 1./ uExact; % finds exact value for R

figure
polarplot(xrange, rEuler, xExactRange, rExact) 
legend('Euler solution', 'Exact solution')
% plots both solutions on same graph

% the exact solution is elliptical  
% the Euler solution is a good approximation because it closely follows
% the exact solution

% following code investigates the effect of changing V0
uExact1 = sym(zeros(1,1000)); 
a1 = subs(A, h, 1); a2 = subs(a1, gam, 2); a3 = subs(a2, V0, 0.5);
% eg. first chooses V0 to be 0.5

for i = 1:1000
    uExact1(i) = subs(a3, x, xExactRange(i));
end
rExact1 = 1./ uExact1;

uExact2 = sym(zeros(1,1000)); 
a1 = subs(A, h, 1); a2 = subs(a1, gam, 2); a3 = subs(a2, V0, 1.3);

for i = 1:1000
    uExact2(i) = subs(a3, x, xExactRange(i));
end
rExact2 = 1./ uExact2;

uExact3 = sym(zeros(1,1000)); 
a1 = subs(A, h, 1); a2 = subs(a1, gam, 2); a3 = subs(a2, V0, 1.45);

for i = 1:1000
    uExact3(i) = subs(a3, x, xExactRange(i));
end
rExact3 = 1./ uExact3;

uExact4 = sym(zeros(1,1000)); 
a1 = subs(A, h, 1); a2 = subs(a1, gam, 2); a3 = subs(a2, V0, 1.6);

for i = 1:1000
    uExact4(i) = subs(a3, x, xExactRange(i));
end
rExact4 = 1./ uExact4;

% finds R for 4 different values of V0

figure
polarplot(xExactRange, rExact1, xExactRange, rExact2, xExactRange, rExact3, xExactRange, rExact4) 
legend('Solution when V0 is 0.5', 'Solution when V0 is 1.3', 'Solution when V0 is 1.45', 'Solution when V0 is 1.6')
% plots R for all different values of V0

% polarplot shows that as V0 increases, the orbits become larger
% so for higher initial radial velocities, the path of the planet around
% the sun will be larger (when 0 <= V0 < sqrt(3))