Looks like we have a 2nd-order differential equation of function . To clean the problem up, let's define some constants:
, , , ,
Then we have the following expression:
Now get by itself so we can use a numerical method to integrate the system:
Suppose you define . Then the system is determined for when you specify the initial conditions: .
Notice that .
Let a very small change in t be denoted . Then you could approximate the value of using the relation .
So the process looks like this:
1. Start at and specify the initial conditions and . This allows you to calculate .
2. Use and to determine and h after a tiny step . Now you can calculate .
3. Use and to determine and h after another tiny step . Update the value of .
4. Repeat until .
This process is called "rectangular integration" and it allows you to numerically solve DEs quickly. If you want higher accuracy, you can use a Runge-Kutta method. These methods follow a very similar procedure but reduce the estimation error at every step. A popular one is a fourth-order method commonly called "RK4."
EDIT: I played around simplifying the original equations. Notice that .