Plotting a Trajectory using Euler's Method

I'm solving the system of ODEs given by

R' = U

U' = Lz^2/R^2 +GM/R^2

using Euler's method. This is no problem, I have a functioning code given by

function euler                                       
Tsim   = 10;           % simulate over Tsim seconds
h      = 0.01;         % step size
N      = Tsim/h;       % number of steps to simulate
x      = zeros(2,N);
t      = zeros(1,N);
G      = 6.67*power(10,-11);
M      = 100;
Lz     = 20004444444; 
R_1    = 10000;
U_1    = 20000;
x(1,1) = U_1;          %IC_ODE_1
x(2,1) = (G*M/(power(R_1,2))+power(Lz,2)/(power(R_1,3)));  %IC_ODE_2

for k=1:N-1
    t(k+1)= t(k) + h; 
    x(:,k+1) = x(:,k) + h* Orbit(t(k),x(:,k));
end

%Compare with ode45
x0 = [1000;0];
[T,Y] = ode45(@Orbit,[0 Tsim],x0);
plot(t,x(1,:)) %Results from Euler algorithm
hold on 
plot(T,Y(:,1),'^') %Results from ode45
legend('Euler','ode45') 
end

function dx= Orbit(t,R)
G  = 6.67*power(10,-11);
M  = 100;
Lz = 200; 
% Equations of motion
dx = zeros(2,1);
dx(1)=R(2);
dx(2)= (G*M/(power(R(1),2))+power(Lz,2)/(power(R(1),3)));
end

This generates plots fine. I want to plot the trajectory however. To do this, I need to be able to generate a value for two new variables. Labelling these by xvar and yvar, I essentially want to use the result of the Euler method code (which gives me a radius at time t) to give me trajectory in polar coordinates. Thus I require something like

theta(:,k+1) = theta(:,k) + h*Angle(t(k),theta(:,k));
    xvar(:,k+1) = x(:,k+1)*cos(theta(:,k+1));
    yvar(:,k+1) = x(:,k+1)*sin(theta(:,k+1));

Where theta is given by the ODE dtheta/dt = Lz^2/R^2. Therefore, I'm assuming I require something like this

function da = Angle(t,R)
Lz = 200;
da = Lz/(power(R(1),2))
end

But theta is given in terms of only one ODE, the one above, but that R(1) is obtained from solving the system above. I just want to be able to plot xvar against yvar with theta varying as it should.

Answer 1

[Swathik Kurella Janardhan]

You can calculate the initial theta(:,k) for an initial value of k by ODE

dtheta/dt = Lz^2/R^2

and then use this value for further calculations of theta(:,k+1) by

theta(:,k+1) = theta(:,k) + h*Angle(t(k),theta(:,k));