Looks like we have a 2nd-order differential equation of function 
. To clean the problem up, let's define some constants:
. 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:
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: 
.
. 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 
.
. 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 
.
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 
.
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 .
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 
.
.
