(para indicar que es la derivada segunda)
Por eso es que podés despejar la acelareción que sufre el choete en la ecuación de la siguiente manera:
y(t)=F/m
o sea
y(t)=(E(t)-P)/m
donde los valores de E(t) son los que obtuvimos en la toma de datos del ensayo del motor.
Teniendo en cuenta que el peso P=mg, en realidad podríamos escribir:
y=E/m - g
Esa es la ecuación diferencial que estoy resolviendo numéricamente. OK, cómo hago eso no lo voy a explicar acá, porque el fundamento es más complicado, pero es una simple receta, una vez que entendiste de donde sale la receta, la aplicás a ojos cerrados y no tiene ningún misterio.
Lo único que determina la solución serán las condiciones iniciales (altura y velocidad inicial) que asumí que ambas valían cero.
Una mejora que puede agregarse es tener en consideración la fricción del aire. Generalmente la fricción es una función de la velocidad: a mayor velocidad, mayor fricción.
Aparecería entonces una nueva fuerza:
F(t)=E(t) - P - f(y')
donde f(y') es la fuerza de fricción, que depende de la velocidad y' (derivada primera de la distancia)
Así la ecuación diferencial quedaría:
y=E(t) - P - f(y')
y se resuelve de manera similar a la anterior. Si la altura que alcanza el cohete es MUY alta, entonces la fricción también dependerá de la altura (ya que varía la densidad del aire y ese es uno de los factores que influye en la fricción, pero para bajas alturas se puede suponer constante y por lo tanto independiente de la altura que alcance el cohete.
Hasta ahi es realmente sencillo de simular porque los datos o los tenes o te los va dando la simulación (como la altura y la velocidad en cada instante). Se me complica un poco la estimación de ritmo de variación de masa del motor, para tenerlo en cuenta de manera de mejorar la simulación, pero ya le voy a encontrar una vuelta. Seguramente se podrá estimar la pérdida de masa en cada instante mediante el cociente de la fuerza de empuje en el instante t dividido la integral de la fuerza de empuje, multiplicado por la masa de combustible. Eso me daría la fracción de masa perdido en ese delta t… suponiendo que la fuerza de empuje es proporcional a la masa eyectada, lo cual no es estrictamente cierto, porque también es propocional a la velocidad de eyección de esa masa, pero más o menos sirve 
Igual no creo que la influencia de la variación de masa sea demasiado importante. Calculo que el término de fricción debe ser más importante y es realmente una pavada de incluir.
Hasta ahi la simulación es 1D… suponemos que el cohete sube derechito y baja derechito (y de culo)… habría que al menos llegar a una simulación 2D, que debe incluir el efecto del viento, o por lo menos un ángulo de despegue.
Con el método de Runge-Kutta, en lugar de resolver la ecuación uno lo que obtiene son aproximaciones numéricas a la solución de la ecuación para ciertos valores de la variable independiente, en este caso el tiempo.
Por supuesto que las condiciones iniciales (altura y velocidad inicial) determinan la solución.
Lo que no aclaré es que al resolver la ecuación diferencial de segundo orden, la solución es la altura y(t). Como el método de Runge-Kutta se puede aplicar solamente a EDOs de primer orden, lo que tenés que hacer para aplicarla a una de segundo orden es aplicar un cambio de variable y transformar la EDO de 2do orden en un sistema de 2 EDOs de 1er orden.
En nuestro caso el cambio de variable lo podríamos escribir así:
y'=v
(lo cual es obvio, pues uno lo que lee es: la derivada primera de la altura es la velocidad)
Por lo tanto y=v' (la aceleración es la derivada primera de la velocidad)
y eso da lugar al sistema de 2 EDOs de 1er orden:
y'=v
v'=E(t)-P
La solución de la primera ecuación es la altura y(t) y le corresponde en nuestro caso la condición inicial y(0)=0
La solución de la segunda ecuación es la velocidad v(t) y le corresponde la condición inicial v(0)=0
Podríamos haber puesto otras condiciones iniciales, pero lo normal en nuestro caso sería que el cohete arranque desde el piso con velocidad cero.
Ahora que tenés eso, le aplicás la “receta Runge-Kutta” que es MUY simple.
Supongamos que tenés una Ecuacion Ddiferencial Ordinaria genérica, expresada de la siguiente manera: y'=f(t,y)
Esto es que el término de la derecha, f(t,y) es una función que depende exclusivamente de la variable independiente (el tiempo) y la solución de la ecuación diferencial, y(t) que es lo que queremos encontrar.
Entonces, si conocés un valor de la solución a tiempo t, aplicando la receta de RK podés encontrar una aproximación a la solución a tiempo t+h, donde h es el incremente, en nuestro caso 1/60 s.
Para eso, hacés el siguiente cálculo (esto tomalo como viene, otro día te explico de donde sale):
K1=h*f(t0,y0)
K2=h*f(t0+h,y0+K1)
y1=y0+(K1+K2)/2
t1=t0+h
¡Listo! t0 e y0 son las condiciones iniciales, con eso encontraste un valor aproximado a la solución para t1=t0+h, y a ese valor le llamás y1.
OK, ¿qué hago ahora?
Pues muy simple, ahora conozco la solución a tiempo t1, puedo volver a aplicar la receta RK para encontrar la solución a tiempo t2=t1+h.
Llamale a eso y2.
Ahora con t2 e y2, vuelvo a aplicar la receta y calculo y3 para el tiempo t3
… y así sigo, todo lo que quiera. Es realmente una papa 8)
Ahora bien, la receta que te pasé es para aplicar a 1 EDO de 1er orden, pero nosotros tenemos un sistema de 2 EDOs de 1er orden.
No pasa nada, supongamos que el sistema de EDOs es:
y'=f(t,y,v)
v'=g(t,y,v)
con condciones iniciales t0, y0 y v0 (fijate que a cada ecuación debe corresponderle una condición inicial de la correspondiente solución).
Le aplicás la receta a CADA ecuación en ORDEN. Para no confundir, identificaremos con K las Kes de la 1er ecuación y con L las Kes de la segunda ecuación. Así que al aplicar la receta queda así:
K1=h*f(t0,y0,v0)
L1=h*g(t0,y0,v0)
K2=h*f(t0+h,y0+K1)
L2=h*g(t0+h,y0+K1,v0+L1)
y1=y0+(K1+K2)/2
v1=v0+(L1+L2)/2
t1=t0+h
¡y listo! Si tenés un sistema de N EDOs… aplicás la receta N veces, una a cada EDO… una verdarera chotada. Claro que si tenés N EDOs lo más probable es que te convenga trabajar matricialmente, porque te podés llegar a hacer un entrevero lindo 
A esta altura, si mirás el programa de AWK que escribí, vas a empezar a reconocer la “receta de RK” por ahi. No pretendas entender al detalle el funcionamiento del programa, porque si sabés algo de programación te va a aprecer que le faltan cosas. Por ejemplo la parte que se lee la tabla de datos debe parecer muy oscura (es una sóla linea), pero eso es porque AWK ya tiene implícito la lectura de datos: los tipos que lo diseñaron sabían bien lo que hacían. AWK es un programa para procesar datos de texto, así que para qué molestarse en que el usuario tenga que programar eso: se da por defecto.
Luego que termina la parte de lectura de datos, el programa está dividido en 2 secciones: mientras tengo datos de la tabla de empuje, calculo las ecuaciones considerando el empuje, cuando se me terminó la tabla de valores de empuje (el motor se apagó) arranca la segunda sección, de vuelo libre. Viste que en el cálculo desapareció el empuje E, y sólo queda el peso P (en realidad como va dividido por m, alcanzaría con poner g es lugar de P/m). Pero viste que el programa es rechoto, igual que la receta.
OK, volviendo al método numérico, RK de orden 2 no es tan bueno como las recetas de RK de orden superior, pero para empezar está bien.
Te paso la receta de orden 4:
K1=h*f(t,y)
K2=h*f(t+h/2,y+K1/2)
K3=h*f(t+/2,y+K2/2)
K4=h*f(t+h,y+K3)
t=t+h
y=y+(K1+2*K2+2*K3+K4)/6
al final soy un nabo, demoro más en escribir toda esta explicación que en reescribir el programa para hacerlo con RK4 e incluir la fricción del aire
Bueno, acabo de hacer la simulación incluyendo la fricción del aire.
La fórmula de fuerza de fricción que estoy usando es:
f(v)=0.5*Cd*da*A*v²
donde
Cd: es el coeficiente de fricción del cuerpo, que en el caso del cohete consideré que valía 0.7, según una publicación de NAR. Este coeficiente es muy variable, depende de varios factores, no sólo de la forma del cuerpo, sino también de la velocidad (a través del número de Mach y el número de Reynolds que mide la turbulencia del fluido que atraviesa), pero para bajas velocidades (subsónicas) se puede considerar independiente de Re y Ma. Acá es donde toma importancia la terminación del cohete: que tan áspera o suave sea la superficie, que elementos extraños perturben el flujo del aire, etc.
da: es la densidad del aire, que tomé como constante y de magnitud 1.29 kg/m³
A: es el área de la sección transversal del cohete, que es Pi*r², medido en m² , por lo tanto tomé r=0.025 m
v: es la velocidad del vector.
Es evidente que esta fuerza varía también en el tiempo, a medida que varía la velocidad del vector.
Así que ahora las ecuaciones quedan así:
y'=v
v'=(E(t)-P-f(v