Ошибка с интегратором Рунге-Кутты 4-го порядка

Я работал над решателем 4-го порядка Рунге-Кутта, и у меня возникли некоторые трудности. Я написал решатель на основе статьи на gafferongames, но когда я запускаю небольшой отдельный пример, ошибки, которые я получаю, намного хуже, чем те, что я получал при простой интеграции Эйлера, даже для простой силы тяжести. Я привел это в отдельный пример (~ 60 строк кода, включая печать), но для его запуска требуется GLM.
Это показывает мою проблему в полном объеме. В строке 55 выводится разница между аналитическим раствором и раствором RK4. Это должно быть относительно мало, но оно взрывается даже после 10 лишних шагов.

#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;
}

Если бы кто-нибудь мог помочь мне здесь, я в конце своей привязи с этим.