因为你是用dvx来计算的，它应该是dvy。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

const int dt = 86400; // 1 day in seconds //
const double G = 6.67e-11;
const double au = 1.496e11;
const double M = 1.99e30;

double vx(double x, double y);
double vy(double x, double y);
double dx(double x, double y, double t);
double dy(double x, double y, double t);
double dvx(double x, double y, double t);
double dvy(double x, double y, double t);

int main(){

  double initial_x = au;
  double initial_y = 0;
  double initial_vx = vx(initial_x, initial_y);
  double initial_vy = vy(initial_x, initial_y);
  double t = 0;

  double x = initial_x;
  double y = initial_y;
  double vx = initial_vx;
  double vy = initial_vy;

  for(int i=0;i<365;i++){
       double k1x = dt* dx(x,y,t);  
       double k2x = dt * dx(x + k1x/2, y + k1x/2, t + dt/2);
       double k3x = dt * dx(x + k2x/2, y + k2x/2, t + dt/2);
       double k4x = dt * dx(x + k3x, y + k3x, t + dt);
       double kx = (1.0/6) * (k1x + 2*k2x + 2*k3x + k4x);

       double k1y = dt * dy(x,y,t);
       double k2y = dt * dy(x + k1y/2, y + k1y/2, t + dt/2);
       double k3y = dt * dy(x + k2y/2, y + k2y/2, t + dt/2);
       double k4y = dt * dy(x + k3y, y + k3y, t + dt);
       double ky = (1.0/6) * (k1y + 2*k2y + 2*k3y + k4y);

       double k1vx = dt * dvx(x,y,t);
       double k2vx = dt * dvx(x+k1vx/2, y+k1vx/2, t + dt/2);
       double k3vx = dt * dvx(x+k2vx/2, y+k2vx/2, t + dt/2);
       double k4vx = dt * dvx(x+k3vx, y+k3vx, t+dt);
       double kvx = (1.0/6) * (k1vx + 2*k2vx + 2*k3vx + k4vx);

       double k1vy = dt * dvy(x,y,t);
       double k2vy = dt * dvy(x+k1vy/2, y+k1vy/2, t + dt/2);
       double k3vy = dt * dvy(x+k2vy/2, y+k2vy/2, t + dt/2);
       double k4vy = dt * dvy(x+k3vy, y+k3vy, t+dt);
       double kvy = (1.0/6) * (k1vy + 2*k2vy + 2*k3vy + k4vy);

       x = x + kx;
       y = y + ky;
       vx = vx + kvx;
       vy = vy + kvy;

       printf("%.3e\t%.3e\t%.3e\t%.3e\n", x, y, vx, vy);
  }

  return 0;
}


// Function for the x velocity of a planet//
double vx(double x, double y)
{
    double theta = atan(y/x);
    double xVel = sqrt((G*M) / (sqrt(x*x + y*y))) * sin(theta);
    return xVel;
}

// Function for the y velocity of a planet //
double vy(double x, double y) 
{
    double theta = atan(y/x);
    double yVel = sqrt((G*M) / (sqrt(x*x + y*y))) * cos(theta);
    return yVel;
}

// Function for dx //
double dx(double x, double y, double t)
{
    double xVel = vx(x,y);
    double dX = xVel*t;
    return dX;
}

// Function for dy //
double dy(double x, double y, double t)
{
    double yVel = vy(x,y);
    double dY = yVel*t;
    return dY;
}

// Function for dvx //
double dvx(double x, double y, double t)
{
    double dVX = ((-G*M*x) / pow(x*x+y*y, 3/2)) * t;
    return dVX;
}

// Function for dvy //
double dvy(double x, double y, double t)
{
    double dVY = (((-G*M*x) / pow(x*x+y*y, 3/2))) * t;
    return dVY;
}

您应该检查的其他一些问题：
要考虑到行星与星星的距离会随时间而变化这一事实
您的代码不会检查行星的位置是否在模拟的范围内

• 你的物理学是错误的。运动方程并不明确地依赖于时间。例如，为什么dvx

double dVX = ((-G*M*x) / pow(x*x+y*y, 3/2)) * t;

在不同的时间应该不同？它应该只取决于x, y。
不要乘以t。

• Runge-Kutta实现也是错误的。

double k1x = dt(x,y,t);  
 double k2x = dt * dx(x + k1x/2, y + k1x/2, t + dt/2);
 double k3x = dt * dx(x + k2x/2, y + k2x/2, t + dt/2);
 double k4x = dt * dx(x + k3x, y + k3x, t + dt);
 double kx = (1/6) * (k1x + 2*k2x + 2*k3x + k);

乘以dt太多。请记住dt是 small。每次乘以一个小值都会降低后续k*x的校正能力。您 * 应该 *

double k1x = dx(x, y, t);
 double k2x = dx(x + dt * k1x/x, y + dt * k1x/2, t + dt/2);
 ....
 double kx = (dt/6) * (k1x + 2*k2x + 2*k3x + k4x);

 // etc

让我们再深入了解一下

double k2x = dt * dx(x + k1x/2, y + k1x/2, t + dt/2);

从变量名可以明显看出y + k1x/2是错误的。

y + k1y/2

这当然需要在计算k2x之前计算k1y。

只有在状态向量的所有分量中完成当前阶段后，才能开始计算Runge-Kutta方法的下一阶段。

顺便说一句正确的方法是

x' = vx
y' = vy
vx' = ax(x,y)
vy' = ay(x,y)

这样k2x的正确行将读取

double k2x = dt * (vx + k1vx/2);

这当然再次要求预先计算k1vx。