Euler's method to plot orbital trajectory of comet

I am trying to use eulers method to plot the trajectory of a comet but no matter what I i get linear plots for position,velocity, and acceleration and generally numbers that don't make sense. I am pretty sure it has something to do with my step size, but Ive tried hundreds f combinations and I dont know how exactly to fix that. The units are all in kilometers.

% Euler's Method
% Initial conditions and setup
h = 1;  % step size
t = 0:h:100;  % the range of x
y = zeros(size(t));  % allocate 
x = zeros(size(t));
vx = zeros(size(t));
vy = zeros(size(t));
ax = zeros(size(t));
ay = zeros(size(t));
y(1) = 0;  % the initial values
x(1) = 193929400;
vx(1) = -5.909;
vy(1) = 50.00294822;
[ax(1),ay(1)] = newaccel(x(1),y(1));

n = numel(y);  
for i=1:n-1
    vx(i+1) = vx(i) + h * ax(i);
    
    vy(i+1) = vy(i) + h * ay(i);
    
    y(i+1) = y(i) + h * vy(i+1);
 
    x(i+1) = x(i) + h * vx(i+1);
 
    [ax(i+1),ay(i+1)] = newaccel(x(i+1),y(i+1));
    
end

This is what my function to determine acceleration looks like.

function [a_newx a_newy] = newaccel(Xx_n,Xy_n)

X_n = [Xx_n * 1000,Xy_n * 1000];
G = 6.67408e-11;
ms = 1.989e30;
mag_X = sqrt(X_n(1).^2 + X_n(2).^2);
a_new = vpa((-(G*ms)/(mag_X)^3)*(X_n));
a_newx = a_new(1)/1000;
a_newy = a_new(2)/1000;
end

[Alan Stevens]

Are your times in seconds (your G value suggests this)? If so, 100 seconds won't take you very far!

[Dominic Troche]

Yes and I've tried higher times, but the loop takes impracticly long (5+ hours) unless I increase the step size to equally as useless scales. What in my code is making anything over 1000 or so iterations take forever?

[Andres Cano]

How did you get Vx? Please answer if you see this

[Andres Cano]

I meant to say how did you get Vy?

Answer 1

[Alan Stevens]

It's the vpa function slowing everything. There is no need for it here. Try the following:

% Euler's Method
% Initial conditions and setup
h = 10;  % step size
t = 0:h:10^7;  % the range of x
y = zeros(size(t));  % allocate 
x = zeros(size(t));
vx = zeros(size(t));
vy = zeros(size(t));
ax = zeros(size(t));
ay = zeros(size(t));
y(1) = 0;  % the initial values
x(1) = 193929400;
vx(1) = -5.909;
vy(1) = 50.00294822;
[ax(1),ay(1)] = newaccel(x(1),y(1));

n = numel(y);  
for i=1:n-1
    vx(i+1) = vx(i) + h * ax(i);
    
    vy(i+1) = vy(i) + h * ay(i);
    
    y(i+1) = y(i) + h * vy(i+1);
 
    x(i+1) = x(i) + h * vx(i+1);
 
    [ax(i+1),ay(i+1)] = newaccel(x(i+1),y(i+1));
    
end

plot(x,y), grid

function [a_newx, a_newy] = newaccel(Xx_n,Xy_n)

X_n = [Xx_n * 1000,Xy_n * 1000];
G = 6.67408e-11;
ms = 1.989e30;
mag_X = sqrt(X_n(1).^2 + X_n(2).^2);
a_new = (-(G*ms)/(mag_X)^3)*(X_n);    % No need for vpa here 
a_newx = a_new(1)/1000;
a_newy = a_new(2)/1000;
end