Error with Runge-Kutta 4th order integrator

Question

I've been working on a 4th Order Runge-Kutta solver, and am having some difficulties. I've written the solver based on the article on gafferongames, but when I run a small contained example, the errors I'm getting are far far worse than what I was getting with simple Euler Integration, even for a simple gravity force. I've tidied it up into a self contained example (~60 lines of code including printing), but it requires GLM to run.
It shows my issue in its entirety. Line 55 prints out the difference between the analytic solution and the RK4 solution. This should be relatively small, but it explodes even after the 10 odd steps it takes.

#include <iostream>
#include <glm/glm.hpp>
struct State{
    glm::vec3 position, velocity;
};
class Particle{
public:
    glm::vec3 position, velocity, force;
    float mass;

    void solve(float dt);
    glm::vec3 acceleration() const {return force/mass;}
    State evalDerivative(float dt, const State& curr);
    void analytic(float t,  glm::vec3 a);
};
int main(int argc, char* argv[]){
    Particle p;
    p.position = glm::vec3(0.f);
    p.mass = 1.0f;
    for(int i = 1; i <= 10; i++)    {
        p.force = glm::vec3(0.f, -9.81f, 0.f);
        p.solve(.016f);
        p.analytic(i*.016f, glm::vec3(0.f, -9.81f, 0.f));
    }
    getchar();
    return 0;
}
void Particle::solve(float dt){
    State t;t.position = glm::vec3(0.f); t.velocity = glm::vec3(0.f);
    State k1 = evalDerivative(0, t);
    State k2 = evalDerivative(dt*.5f, k1);
    State k3 = evalDerivative(dt*.5f, k2);
    State k4 = evalDerivative(dt, k3);

    position += (k1.position + 2.f*(k2.position + k3.position) + k4.position)/6.f;
    velocity += (k1.velocity + 2.f*(k2.velocity + k3.velocity) + k4.velocity)/6.f;
    force = glm::vec3(0.f);
}
State Particle::evalDerivative(float dt, const State& curr){
    State s;
    s.position = position + curr.position*dt;
    s.velocity = velocity + curr.velocity*dt;

    s.position = s.velocity;
    s.velocity = acceleration();
    return s;
}
void Particle::analytic(float t, glm::vec3 a){
    glm::vec3 tPos = glm::vec3(0.f) + 0.5f*a*t*t;
    glm::vec3 tVel = glm::vec3(0.f) + a*t;

    glm::vec3 posdiff = tPos - position;
    glm::vec3 veldiff = tVel - velocity;
    std::cout << "POSITION: " << posdiff.x << ' ' << posdiff.y << ' ' << posdiff.z << std::endl;
    std::cout << "VELOCITY: " << veldiff.x << ' ' << veldiff.y << ' ' << veldiff.z << std::endl << std::endl;
}

If anyone could help me out here, I'm at the end of my tether with this.

Answer 1

Well, I feel stupid. I've been working on this for a few hours now, and I missed a little step:

position += (k1.position + 2.f*(k2.position + k3.position) + k4.position)/6.f;
velocity += (k1.velocity + 2.f*(k2.velocity + k3.velocity) + k4.velocity)/6.f;

should be

position += (k1.position + 2.f*(k2.position + k3.position) + k4.position)*dt/6.f;
velocity += (k1.velocity + 2.f*(k2.velocity + k3.velocity) + k4.velocity)*dt/6.f;