Unity freezes when I run this and I have no idea why, please help

Hey everyone, so I am trying to do an orbit simulation with multiple bodies involved, currently there are nine and I let the simulation run for a while and then it completely freezes and I have to close the program in order to work on it. I am not sure how to fix this, the code is pasted below.

using System.Collections; using System.Collections.Generic; using UnityEngine; using System;

public class PlanetPropagator : MonoBehaviour {

 public float semiMajorAxis;
 public float inclination;
 public float eccentricity;
 public float rightAscension;
 public float argOfPerigee;
 public float meanAnomaly;
 public float meanMotion;
 float t;
 public KeplerianMotion keplerianMotion;
 float[] planetPosition = new float[3];

 // Use this for initialization
 void Start () {
     DateTime satelliteEpoch = new DateTime(2015, 4, 20, 22, 45, 0);
     DateTime tVernal = new DateTime(2015, 3, 20, 22, 45, 0);
     t = (float)(satelliteEpoch - tVernal).TotalSeconds;

     keplerianMotion = new KeplerianMotion();
     planetPosition = keplerianMotion.GetNewPosition(eccentricity,meanAnomaly,meanMotion,semiMajorAxis,inclination,rightAscension,argOfPerigee,t);
     transform.position = new Vector3(planetPosition[0], planetPosition[1], planetPosition[2]);
     t += 1000;
 }
 
 // Update is called once per frame
 void Update () {
     planetPosition = keplerianMotion.GetNewPosition(eccentricity, meanAnomaly, meanMotion, semiMajorAxis, inclination, rightAscension, argOfPerigee, t);
     Debug.Log(planetPosition[0] + ", " + planetPosition[1] + ", " + planetPosition[2]);
     transform.position = new Vector3(planetPosition[0], planetPosition[1], planetPosition[2]);
     t += 1000;
     //Debug.Log("Time: " + t);
 }

}

public class KeplerianMotion { private static float muEarth = 398600.4418e+9f; private static int secsPerDay = (24 60 60);
private static float PI = 3.14159265359f; private static float deg2rad = PI / 180f; private static float wEarth = 7.29211e-5f;

 public float[] GetNewPosition(float eccentricity, float meanAnomaly, float meanMotion, float semiMajorAxis, float inclination, float rightAscension, float argOfPerigee, float t)
 {
     float[] orbElements = new float[7];
     float a, e, i, ra, arg, ma, mm;
     float eOld;
     float r;
     float TA;
     float F_e, dF_e, eNew;
     float[] pos = new float[3];
     float[] posECI = new float[3];
     float[] tempPos = new float[3];
     float[,] matrix = new float[3, 3];
     float[] posArray = new float[3];

     //DateTime satelliteEpoch = new DateTime(2015,4,20,22,45,0);

     e = eccentricity;
     ma = meanAnomaly * deg2rad;
     mm = meanMotion * deg2rad;
     a = semiMajorAxis;
     i = inclination * deg2rad;
     ra = rightAscension * deg2rad;
     arg = argOfPerigee * deg2rad;

     //DateTime tVernal = new DateTime(2015, 3, 20, 22, 45, 0);
     //float t = (float)(satelliteEpoch - tVernal).TotalSeconds;

     float n = mm * 2 * PI / secsPerDay;
     meanAnomaly = meanAnomaly + n * t;
     eOld = meanAnomaly;

     // Newton Raphson formulation
     F_e = eOld - e * Mathf.Sin(eOld) - meanAnomaly;

     dF_e = 1 - e * Mathf.Cos(eOld);

     eNew = eOld - F_e / dF_e;

     //Iteration until convergence
     while (Mathf.Abs(eNew - eOld) > 1e-5f)
     {
         eOld = eNew;

         F_e = eOld - e * Mathf.Sin(eOld) - meanAnomaly;

         dF_e = 1 - e * Mathf.Cos(eOld);

         eNew = eOld - F_e / dF_e;
     }

     // Returning the out values
     TA = 2 * Mathf.Atan(Mathf.Sqrt((1 + e) / (1 - e)) * Mathf.Tan(eNew / 2));

     //Debug.Log("TA: " + TA);

     r = a * (1 - Mathf.Pow(e, 2)) / (1 + e * Mathf.Cos(TA));
     //Debug.Log("R: " + r);
     posArray[0] = r * Mathf.Cos(TA);
     posArray[1] = r * Mathf.Sin(TA);
     posArray[2] = 0;

     //Debug.Log("posArray" + posArray[0] + ", " + posArray[1] + ", " + posArray[2]);
     float SR = Mathf.Sin(ra), SA = Mathf.Sin(arg), SI = Mathf.Sin(i);
     float CR = Mathf.Cos(ra), CA = Mathf.Cos(arg), CI = Mathf.Cos(i);
     //Debug.Log(SR + ", " + SA + ", " + SI + ", " + CR + ", " + CA + ", " + CI);

     matrix[0, 0] = CR * CA - SR * SA * CI;
     matrix[0, 1] = -CR * SA - SR * CA * CI;
     matrix[0, 2] = SR * SI;
     matrix[1, 0] = SR * CA + CR * SA * CI;
     matrix[1, 1] = -SR * SA + CR * CA * CI;
     matrix[1, 2] = -CR * SI;
     matrix[2, 0] = SA * SI;
     matrix[2, 1] = CA * SI;
     matrix[2, 2] = CI;

     for (int ii = 0; ii < 3; ii++)
     {
         for (int j = 0; j < 3; j++)
         {
             tempPos[ii] = tempPos[ii] + (posArray[j] * matrix[ii, j]);
             //Debug.Log("tempPos: " + tempPos[ii]);
         }
     }
     posECI[0] = tempPos[0];
     posECI[1] = tempPos[2];
     posECI[2] = tempPos[1];

     return posECI;
     /*
     float[] posECEF = new float[3];
     float ang = wEarth * t;
     float x = posECI[0], y = posECI[1], z = posECI[2];

     posECEF[0] = x * Mathf.Cos(ang) + y * Mathf.Sin(ang);
     posECEF[1] = z;
     posECEF[2] = -x * Mathf.Sin(ang) + y * Mathf.Cos(ang);

     return posECEF;*/
 }

}