Numerical second order pde solver

[DivergentSpectrum]

Edit:whoops wrong forum mods please move

2nd edit: I just had dinner then got back on the computer, input some points and saw a beautiful elipse.(complete with a fascinating flower petal design due to inaccuracies) Weird lol! No idea why it wasnt working before Now to implement RK4 bwahahaha

im writing a program in WPF(an amalgam of c# and xaml) that tracks the position of a particle after a certain amount of time.
the particle is placed in a force field, where force is a function of position and velocity.
I used a crude numerical method, which seemed pretty obvious to me but for some reason it won't work. the particle appears to travel in a straight line as if there was no force.

the number of increments and the final time to be calculated for is decided by the user, and the variable dt is calculated from that.
vx,vy,vz are the velocity vector
ax,ay,az are functions, and also the acceleration vector.

heres the relevant code. I am pretty sure the problem is in the first block.

      while (i < n - 1)
        {
    

            vx[i+1] = vx[i] + dt * ax(x[i], y[i], z[i], vx[i], vy[i], vz[i]);
            vy[i + 1] = vy[i] + dt * ay(x[i], y[i], z[i], vx[i], vy[i], vz[i]);
            vz[i + 1] = vz[i] + dt * az(x[i], y[i], z[i], vx[i], vy[i], vz[i]);


            x[i + 1] = x[i] + vx[i + 1] * dt;
            y[i + 1] = y[i] + vy[i + 1] * dt;
            z[i + 1] = z[i] + vz[i + 1] * dt;

            

            i++;
        }

for verifiability, I am using a gravitational field, also for simplicity i have eliminated force from the equation and I am working with acceleration alone(ill add that later)
here i have methods define the acceleration vectors:

     //acceleration vectors x componentvvvvvvv
        static double ax(double x, double y, double z, double vx, double vy, double vz)
        {
            double ax =  10000*(-x / Math.Pow(Math.Sqrt(x * x + y * y + z * z), 3));
            return ax;
        }
        //acceleration vectors x component^^^^^^


        //acceleration vectors y componentvvvvvvv
        static double ay(double x, double y, double z, double vx, double vy, double vz)
        {
            double ay =  10000*(-y / Math.Pow(Math.Sqrt(x * x + y * y + z * z), 3));
            return ay;
        }
        //acceleration vectors y component^^^^^^

        //acceleration vectors z componentvvvvvvv
        static double az(double x, double y, double z, double vx, double vy, double vz)
        {
            double az =  10000*(-z / Math.Pow(Math.Sqrt(x * x + y * y + z * z), 3));
            return az;

        }
   //acceleration vectors z component^^^^^^

 

[rcgldr]

If you don't want to do RK4, you could try one trapezoidal corrector step:

    vx[i+1] = vx[i] + dt * ax(x[i], y[i], z[i], vx[i], vy[i], vz[i]);
    /* ... */
     x[i+1] =  x[i] + dt * vx[i];
    /* ... */
    vx[i+1] = vx[i] + .5 * dt * (ax(x[i], y[i], z[i], vx[i], vy[i], vz[i]) +
              ax(x[i+1], y[i+1], z[i+1], vx[i+1], vy[i+1], vz[i+1]));
     x[i+1] =  x[i] + .5 *dt * (vx[i]  + vx[i+1]);

update - see that you plan on using rk4, so no need to implement the corrector step method.

 

[AlephZero]

The best way to find the problem this would be to print out all the variables for one time step (or inspect them with your IDE's debugger, if you can believe what it's telling you!) and check the values with a calculator.

The code looks sensible, but if this is a snippet from a system that is built with tools that many people here have probably never used, the cause could be almost anything.