Bas, you have a sign error in the gravity term. It should be:
Use one of matlab's 'ode' solvers to numerically solve these equations. Read about how to use them at:
Note that you would have numerical accuracy difficulty if the angle phi approaches near zero because of the division by sin(phi) in the first equation. At that point theta would change very rapidly. That is inherent in the physical situation as measured by the two angles.