{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# N-body simulation\n",
    "\n",
    "make a simple $N$-body simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "######## EULER #########\n",
    "def integrator_euler_part1(n_particles, x, v, timestep):\n",
    "    for i in range(n_particles):\n",
    "        x[i][0]  += timestep * v[i][0]\n",
    "        x[i][1]  += timestep * v[i][1]\n",
    "        x[i][2]  += timestep * v[i][2]\n",
    "        \n",
    "def integrator_euler_part2(n_particles, v, a, timestep):\n",
    "    for i in range(n_particles):\n",
    "        v[i][0] += timestep * a[i][0]\n",
    "        v[i][1] += timestep * a[i][1]\n",
    "        v[i][2] += timestep * a[i][2]\n",
    "\n",
    "        \n",
    "        \n",
    "######## LEAPFROG #########\n",
    "def integrator_leapfrog_part1(n_particles, x, v, half_timestep):\n",
    "    for i in range(n_particles):\n",
    "        x[i][0] += half_timestep * v[i][0]\n",
    "        x[i][1] += half_timestep * v[i][1]\n",
    "        x[i][2] += half_timestep * v[i][2]\n",
    "\n",
    "def integrator_leapfrog_part2(n_particles, x, v, a, timestep, half_timestep):\n",
    "    for i in range(n_particles):\n",
    "        v[i][0] += timestep * a[i][0]\n",
    "        v[i][1] += timestep * a[i][1]\n",
    "        v[i][2] += timestep * a[i][2]\n",
    "        x[i][0]  += half_timestep * v[i][0]\n",
    "        x[i][1]  += half_timestep * v[i][1]\n",
    "        x[i][2]  += half_timestep * v[i][2]\n",
    "\n",
    "        \n",
    "        \n",
    "######### Calculate Gravity, direct method ########\n",
    "def gravity_calculate_acceleration(n_particles, m, x, a, G):\n",
    "    for i in range(n_particles):\n",
    "        a[i][0] = 0.0 # reset\n",
    "        a[i][1] = 0.0\n",
    "        a[i][2] = 0.0\n",
    "        for j in range(n_particles):\n",
    "            if (j == i):\n",
    "                continue\n",
    "            else:\n",
    "                dx = x[i][0] - x[j][0]\n",
    "                dy = x[i][1] - x[j][1]\n",
    "                dz = x[i][2] - x[j][2]\n",
    "                r = math.sqrt(dx*dx + dy*dy + dz*dz)\n",
    "                prefact = -G*m[j]/(r*r*r)\n",
    "                a[i][0] += prefact * dx # sum for all j\n",
    "                a[i][1] += prefact * dy\n",
    "                a[i][2] += prefact * dz\n",
    "                \n",
    "# main        \n",
    "def integrate(m, x0, v0, tmax, dt, G=1.0):\n",
    "    npar = len(m) # number of particles\n",
    "    x = x0\n",
    "    v = v0\n",
    "    # \"create\" an empty list with size npar x 3 (ax,ay,az)\n",
    "    a = [ [0.0 for _ in range(3) ] for _ in range(npar)] # stupid python\n",
    "    \n",
    "    tstep = dt\n",
    "    half_tstep = 0.5*tstep\n",
    "    t = 0.0\n",
    "    xsave = []\n",
    "    \n",
    "    # main iteration\n",
    "    ## backward halfstep for leapfrog\n",
    "    # gravity_calculate_acceleration(npar, m, x, a, G)\n",
    "    # for i in range(npar):\n",
    "    #     v[i][0] -= half_tstep * a[i][0]\n",
    "    #     v[i][1] -= half_tstep * a[i][1]\n",
    "    #     v[i][2] -= half_tstep * a[i][2]\n",
    "    \n",
    "    while (t < tmax):\n",
    "        #integrator_leapfrog_part1(npar, x, v, half_tstep)\n",
    "        integrator_euler_part1(npar, x, v, tstep)\n",
    "        gravity_calculate_acceleration(npar, m, x, a, G)\n",
    "        #integrator_leapfrog_part2(npar, x, v, a, tstep, half_tstep)\n",
    "        integrator_euler_part2(npar, v, a, tstep)\n",
    "        t += tstep\n",
    "        \n",
    "        # save position \n",
    "        # if we append list directly, the value is mutable... f*ck python\n",
    "        # same value for all -_-!\n",
    "        xsave.append([])\n",
    "        for particle in x:\n",
    "            xsave[-1].append([])\n",
    "            for xyz in particle:\n",
    "                xsave[-1][-1].append(xyz)\n",
    "\n",
    "    return xsave\n",
    "\n",
    "                \n",
    "def plot_xy(pos):\n",
    "    '''function to plot xy-coord'''\n",
    "    \n",
    "    npar = len(pos[0])\n",
    "    pos = np.array(pos) # python-list to numpy-array\n",
    "    # for easy slicing\n",
    "    \n",
    "    fig = plt.figure(figsize=(5,5))\n",
    "    ax = plt.subplot(111)\n",
    "    for i in range(npar): \n",
    "        # plot each object with different color\n",
    "        plt.plot(pos[:, i, 0], pos[:, i, 1], ',')\n",
    "        \n",
    "    plt.axis('scaled') # make it using same scale for x and y\n",
    "    plt.xlabel(\"x\"); plt.ylabel(\"y\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test for two-body problem\n",
    "\n",
    "there is analytical solution for this problem.\n",
    "\n",
    "**initial condition**:\n",
    "- $m_1 = 0.6$; $m_2 = 0.4$\n",
    "- $x_1 = (-0.4, 0, 0)$; $x_2 = (0.6, 0, 0)$\n",
    "- $v_1 = (0, -0.2, 0)$; $v_2 = (0, 0.3, 0)$\n",
    "\n",
    "G = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**dt = 0.01, tmax = 100**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAADVCAYAAAAM/NhBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX3MbVdd579fW0iZQdpbei2X3uKDpKjXIi8+EtRJRCih\nL8aLL39IHe0kmKaOKMbXazDmmhmGimZCTBiaOzCxRCtBJbaBa0hbUZPhxT5Xa2mp5VamDhdb+lix\nqEnVys8/nn2O++yz3vbrWvuc7yc5Ofvlt9f6rX3O/p7f+q2196GZQQghxHB8RW4HhBBi05CwCiHE\nwEhYhRBiYCSsQggxMBJWIYQYGAmrEEIMjIRVCCEGRsIqhBADI2EVQoiBOT+3A2NwySWX2M7OTm43\nhBAbxpkzZ/7GzA7H7DZSWHd2drC3t5fbDSHEhkHyr1LslAoQQoiBkbAKIcTASFiFEGJgJKxCCDEw\nElYhhBiYjZwVMHd2Tnx4sLIeufm6wcoSQqQhYc3MkCLapnwJrhDjoVRABnZOfHj5yulDKb5sHScv\nDL/HbGLH18sRWVDEOiGli1fdP0W0HTl5IXDyyX9/X2xbsXly1d71HtuXUn+9PpdfYjQUsU5E6aLa\nRFFsAk0RbIrXSnT55KrQOre3eHfW1RDMtfWG6CrCHQ1FrCNTqjg9cvN1Sb41bbY2knVFoL7IsG7v\njESb5bi6+wHxdL2nHudr11pZimr7kDViJXk1yYdIPkzyhGP/D5C8j+SnSH6M5Etz+NmFkiO+VFF1\niWjJ7RqFtUi0IaquKNK33ixzuRyJSBd29feQv2v+N+v2RbguW0W1XcgmrCTPA/AuANcAOAbgjSSP\nNcz+H4BvN7OXAPhvAE5N62V7ShaehVDG/Euxqw9+bRxrYhbJWS5sU6LH5fEJ3Xqf6LrKb9r4BNXX\nzpDIhqJv4YRmlqdi8lsAnDSz11frPw8AZvZ2j/0hAPeb2WWxsnd3dy3H063mLjKpkazv2NkSEopm\n9zi2vFhvluvritf3tbEZ+j1UR0rbtyR1QPKMme3G7HKmAi4D8Lna+rlqm483Afj9UT3qSMlRahv6\ntGF258AXrXmjOs9IvzcF0Bic8qUF6j60iVBj+dNQCsBXp+scOX8kPNGsItkls5gVQPI7cCCsPxew\nuZHkHsm9/f39yXyblZhMQPECu5YzTYy+fIM7IaFaE9mQgMaEMpJ3bZbT9LtZlvfchM5BQlSqdAGA\nvML6eQCX19aPVttWIPmNAN4D4LiZPeErzMxOmdmume0ePhx9wPcglCggpXTJixNY38BTbNvy+MSu\nvysqbPrhsm1uSxXutXb2SAlE64v86DjbsZ0Cm1NY7wFwBckXknwmgO8HcEfdgOQLAHwQwA+a2Wcy\n+OilKNGoUZJfi5xtNp+W0akjilzsT9rmEMDF9uAxrgjWM8DkFeeA30t/QrnRFuuxKHrt3CbmVbdQ\nYLMJq5k9DeDNAD4C4EEAHzCzB0jeRPKmyuwXATwXwP8ieS/JIv5vpQTxGioyDZXTdd+C+nma7Jw5\nu8eJouQsJyGCddl5UwCebrwvHeFrz1q5ngi02VbX+to5cYhgLIJNFdotEdisOVYzO21mLzazF5nZ\n26ptt5jZLdXyD5vZITN7WfWKjsaNTQmiCqT50Vb8hto35DGtCEVtK3YOwW3ua+4PLdfrCKUJQnXE\nhDTUdW+Tm/WWE4ham7j8Ddk0t21BBDuLwSvRjVJ+BOqMkhpwRYJrNompgLp9iqiGotamP75Ir7nN\n5Wcb4ayXF0sRLMtJTA248EWtsSh2gwVWwtqCEoWqNFJTFL3PZbCL7InGVo5z7Rso17rma6Tr7MzH\nBrr+rvbHBHfFtl5eQEBjIpuSEvAKeNNmswRWwjpzhsi1DjmToI1g9hLXWDTmO6Z+rKu8ul29LN9+\nbx2+9IMvqnT9UDTb5YqIEyPVFV88KYDWIhtY9tUXYoPEVcKaSKnR6hB+5X4ubHL9MeEE/ILhOi4l\nQg36kZBfXRNOT061aV+vKyrKDvHtmm9tI7Kuc1qvK4TvHG9I9CphFUmkRLV9It+ouKYIwcLOeXwk\nWm3ui0WKrvLXRPJC97amny4BXvErVVx9whiJYqPnzuFTs+0hP13rze0bJrAS1gRKjVZDDH2jwOIc\nhMrte56c0Wubi6tvtOorMya2oWNW/GqKTSPX6bP3imdTLH3C6Tkfzqg0IWptkxJIFldXyiJwjgtH\nwrqh+ISwr+BO/iOTEp3GLr5Qt9QXgbUVdF9Xf802JDxPesryiWckDeCN7h1i1yo1EBFXX3tDdot9\ni3asrM9PXPWg6whziFZjYtlWXHdOfLjXk676sqg/fBEGBNdHmyizua31ckPgvD458pVtBDUktE0/\nnHW0SBH4yve1qasghsqcCYpYZ8QjN1/nfI1Rj6++KXjkgut7XJQdBNbZPa93jROWXWWt1RnoIq/Y\nOoSxy3K9rSmpgeXxnhRB00dfmsC13LRzlenDlzIpmGzPYx2TIZ/HmjNiLeWBKjGGPEePXHA9dp66\nbeV9hdRua9M2N6HIdGnzZFyEQpGy7/g2Xfk276G6Y8uuNrch02ea+jxWpQICTC2qcxHSJnW/Oz8o\nuxJSJ10uPl+klQtfVNuMxkJpgeZ6LNfaXF7W0zI1EHz3pD58yynnKcmukB9MD0oFZGbqbvbY9Ekd\nrEWnFenzXAsT0zYsUgW+iHNpF0kB+IR0xSaSGghF+k0h9fkaEkdnHjfhOJ8vBSJhzcQmiWmIer7W\naxOKVkPMWUh9rAmeZ1DJK64h0XQJaGK+NSSkPkFN+XEIbgt8rl2EeEKUY/UwRhpgG4Q0Rv28+vKp\noW0bJaJd8XX5Q/tS8pxd8q1O/wbIsabU07SfgDn859XWsC3RaQptzkUzkt156rZZTH8bnUWkubLu\nyNG60gvRNEDLPOuiDFd5dT+ay9E2OlIOvnNRty8ECeuISFD9NPOpvvwqgLU0gcS1Rl3cnCIaikRd\naYCug1iBrn1IXFPE1pVecO4rR1yVCnAwxIUrQfUQGMFudv1D6PxGcJ3nxfYVu0AaIFh+oKuekgpI\n7eKn4mrrCMwiFUDyapIPkXyY5AnH/q8j+XGS/0Typ3P42BZFqS2pXVzLrn/ChaGoNUIsVVDfDnjS\nAM3Bqicb5UVSC8tyPbMLfJ9zF2EcS7A7kk1YSZ4H4F0ArgFwDMAbSR5rmP0tgB8H8KsTu9caCWoC\nCVFQm1tpJa4JeCPSxHzryjEdUwH1bYt6fDYuu+ZysL35RRXImAog+S0ATprZ66v1nwcAM3u7w/Yk\ngH8wsySB7ZMK6HKxSlATiI0GNy6c1M9B574DvsGsWL41Wm7ELpQKGKqOpu3AzCEVcBmAz9XWz1Xb\nOkHyRpJ7JPf29/d7O5eKLuyOrFxg6xdASg8g54NiZo1r4Co2sBVKCyzL9fxQpghcamrANzXLV2Ym\nNmZWgJmdMrNdM9s9fPjw6PWp69+C0HzHCCnPf5W4dsQlpIBHaGv7fGW5yl2xaUTHa8dFUgPedgQG\nrTKJa05h/TyAy2vrR6ttxSNBbUEPUV2g8z0izilRgRzscn+LPKtLUJs29fVmfcltGX9WQCo5c6zn\nA/gMgNfiQFDvAXC9mT3gsD2JCXKsKZGPLvIWDCCqdWKfjz6bAYjlPX2iGjrGW1fC9Ksh8rIL2wEo\nPsdqZk8DeDOAjwB4EMAHzOwBkjeRvAkASD6P5DkAPwngF0ieI/mcXD7rwu2JJ5+ais7/BCyjVdcc\n2Eik6svHNstYbgtEsnUbl33omFg5E5A1x2pmp83sxWb2IjN7W7XtFjO7pVp+zMyOmtlzzOyiavlL\nOXzVRV0G+hwmxDVbICXv6ksB+FIGvnqd+8qYThVjYwavRGG4upAD4RNXDWKNgGuGQHP/ctnVlU8c\nxGrur5edNKvANXCVL9cqYU1EF20PRviCS1wnxjcI5ZuK1ezCr9mmzC5oRMRrAp0g6qH1EZGwiuHp\nM7LbAqUFJqbrDAGX4PrmuQbns+YTyrZIWCtSIh1FQwlM/GV3ias+p5GJ5VtT57G6ltf2Jf4wJ9tN\n8/2UsIrxyJTjUiQ7AbGZAUB4ECu0vFh3bk+NZvPOZZWwtkTRUCITfrGbQqrPaCJSpj758qvRZd9N\nA4lpplBkOkHUKmEVw1FwzkuMSH1+cijC9M1JjQ2KeeudJpffBQlrBxQRRSjsSy4mwncrbH3/yrpr\n5oArHZAy3aqsH3UJq9gIlA4oCG/+NTAo5Yt4F9sG93FcIZawdkQXboMCHoChQauCiIlrcCZA4nep\n4J6RhLUHElchAvQRV+f6BJHsQEhYxXAUHEGITIRuhwXWo9PQg1YKFtImElb0izwVtZaF0gEF4otO\nV2w6pAIKRsIq+jOjSEJkInSr63I9Jr6heasdHsIy4vdWwjoAilox6+hCTERolsDSJjA3NXXSf+xB\nLxMgYR0IiasQCcTEtes81pU68vegsgoryatJPkTyYZInHPtJ8teq/feRfEUOP1ORuAqRQOz5AQub\npX2CwPqeppWJbMJK8jwA7wJwDYBjAN5I8ljD7BoAV1SvGwG8e1InO7B14lpAdCBmSDRybXmjgO/5\nr5nIGbG+EsDDZvZZM/tnAO8HcLxhcxzA++yATwC4iOSRqR1ty9aJqxBdcInrYrtrfUbkFNbLAHyu\ntn6u2tbWpkgkrkIkkHqjQN3eV07MZkI2ZvCK5I0k90ju7e/v53YHgMRViCRc/yzg/TuWyOMJQzYT\nklNYPw/g8tr60WpbWxsAgJmdMrNdM9s9fPjwoI72QeKaB90oMDOiD7z2PLSlUHIK6z0AriD5QpLP\nBPD9AO5o2NwB4Ieq2QGvAvCkmT06taN9kbgKkUD9sYOL9eW+yP9dFSa02YTVzJ4G8GYAHwHwIIAP\nmNkDJG8ieVNldhrAZwE8DOB/A/ivY/gyRXQjcRWiJc5nCyT+NUtmsuZYzey0mb3YzF5kZm+rtt1i\nZrdUy2ZmP1rtf4mZ7eX0ty8S1/HROd4A+j7RqoDodWMGr+bCzokP6+IfGeVXN4A+U64KiF4lrGJj\n0A/WhpHyZ4WFImHNxMZFrjP74osZUUDXvi0S1opc3cdZi2tBYro4j0oDbBiuB2TPQGglrAUw2+i1\nkC/4LM+dSCc2gFXI97COhLXGmNFOStmzFdgFBUWwYkNxPtC64/duREGWsE5EG8GctbhOTP1cKQ2w\nwRTw8Oo2SFgblHJxzj56nQCJ6pZRYJffR1RYSf4YyUNTOCPWKV5gM33Ziz4nYutJiVgvBXAPyQ9U\nT/zn2E7lpsTop3iBnZDmeSjx8xIjMZOoNSqsZvYLOHiC/3sB/BcAZ0n+D5IvGtk30eCRm68rT1xj\nD8cYGImqmANJOVYzMwCPVa+nARwC8Dsk3zGib1kp8YJdiEpR0etEEURRbRZ5mUHUmpJjfQvJMwDe\nAeD/AniJmf0IgG8C8L0j+5eVEsV1wSJ6LVJwBo5afe0r+fMRI5PyTwIZSYlYLwbwPWb2ejP7bTP7\nFwAwsy8D+M5RvSuAUi/epthkFdgRv8zF/WiIsilkOhYPevmbxe7uru3tDfuEwTle4JP9KHj/LqOb\n4IYiVN26KlboI6Qdvp8kz5jZbtROwprOHMV1wehC1OahxA7anFuJqlgy8V1XqcKa5QYBkheTvJPk\n2erdOU+W5P8h+TjJ+6f20UVJF3RbX5r52FF+JDrcblj3KaVNJX0GogAKyak2yXXn1QkAd5vZFQDu\nrtZd/DqAq6dyak7EhDEkQM3ZBYPkZ11/Y+yot4+4S1TFXMiSCiD5EIBXm9mjJI8A+EMz+1qP7Q6A\nD5nZlanlj5UKqFNCWmCqea2LeqLvF1yPnaduwyMXXA8Aa8t96hfCS5d0wMipgFzC+ndmdlG1TABf\nXKw7bHdQoLAC04hrkTcFBPCJa0xYfe2UqIok2orrXHOsJO8ieb/jdbxuV9180FvdSd5Ico/k3v7+\nft/iktC/u66yENPFu2ufD4mq2CTOH6tgM7vKt4/kF0geqaUCHh+gvlMATgEHEWvf8lJZXPxzEsCx\nWESl9ag1JqguJKhiVCYY8Mo1eHUHgBuq5RsA3J7Jj8GYWgxKFR9X1FpPA4SE9pGbryu2XaJwCpsd\nkEtYbwbwOpJnAVxVrYPk80meXhiR/C0AHwfwtSTPkXxTFm8TmVIYSo2Qd566bSmkzbxqM/+63C5B\nFRtGFmE1syfM7LVmdoWZXWVmf1tt/2szu7Zm90YzO2JmzzCzo2b23hz+tmXbhaIpoL7Bq20+R2Jg\nCrmVdYH+QWBEhhDYOYpPyjSrpejOsH1CxJCwTkAfgS21y59CfarVShRbz4cVFmmImXHywnbfoYly\nsaPNChDr1MV1zoKZQjMVAFTtP5nRKSEmQhFrJsbIwxbZrV7+u2btltdFhNF8FyKVtpHqxEhYM7MQ\n2FRRTHkGQA4Wfq20xyWmC+qCu9xW7oUiCiLyTIoS0GMDCyYklKXc6pr0g7CMTBtC6xLchZ0QLvr+\n+Pb8bqXe0qoca8E0RasupDlEtVeqwSei9X3NNIEEVtQJ9X5SBHfC75OEdUb4hK3+hKmmfdsHSC/K\nGpS1SNUTtSoVIFzUvxfOHk553xulAsS0xKIO17vYXkKimRqp1u17kv3pVkKssRahNmYMLGxcx4nt\nI/a5F/y9kLCK6ViJRl2pAM/UrLaRiZg3vqlUfSLOiXs+ElYxLa4868r2yDxXCexmE0r/zOizl7CK\nPDSj0FD0WrfXYNdmUo9SCxvh74IGr0QeUqbMpKYCCr/IRIAhB6dC5QyEBq9E2azcceWYAaC7tjaf\n1MGpGf5wSlhFPkKDWG2nZQFKEcyF1k+kqqeM5iGyElaRl6RBrBZ5Vs0iKJchHpwy0b+x9iWLsJK8\nmOSdJM9W74ccNpeT/CjJT5N8gORbcvgqRqaZBogJbRvRLPwJSFtB7DNwpXU2gFwR6wkAd5vZFQDu\nrtabPA3gp8zsGIBXAfhRkscm9FFMhWuGQH3ZKa6OWQP1slz5V4nsdKyM8oeeE9FM43g+15mRS1iP\nA7i1Wr4VwBuaBmb2qJn9abX89wAeBHDZZB6KaWkOXC3evbe6+gTXczts81gxPEP8gMXy5W2ENqMo\n5xLWS83s0Wr5MQCXhoxJ7gB4OYBPBmxuJLlHcm9/f38oP8WUuNIBoZxqG3F1pRya5YluhH6shs53\nz+TzGu3pViTvAvA8x6631lfMzEh6J9OSfDaA3wXwE2b2JZ+dmZ0CcAo4mMfayWmRn6a41pd9zxFI\nEVdXGbEoV7iZSkRnzGgRa/W31lc6XrcD+ALJIwBQvT/uKoPkM3Agqr9pZh8cy1dRGPVocmW5kSJY\nLrcUV1+dShX4cZ2X5rlMuemjlw/zSAMA+VIBdwC4oVq+AcDtTQOSBPBeAA+a2f+c0DeRm2VuNXDD\nQFNwu4prcNDswvW6tgHXDxsQFso1ke0hqq58e92vGZBLWG8G8DqSZwFcVa2D5PNJnq5svg3ADwJ4\nDcl7q9e1edwVWfClA1yRUpd0gPMOLs82l+gulueO74cK8AvkkNHo2kwAz4/djNCzAsQ8aNPNjOZW\nE25IiAlMbN/C51JytqHzAcS3pZ6L0LY1n0bKyY54zvWsALE5rHXfXSP+vqlVgalXLtvmNm99LntH\nxLXc74ly+0a+rnJd7z7xWytvwh+CNm0u5QcqEQmrKB9XV9Gbc3VMp3J2/xMjVmcZPntXdOcQ5Hr9\na21NEMnm8cF0RygK9fi7tm3gaDVVJLvkVwsRYAmrmA9NsWtuq6/HItQ1+4SIdXmMw967zxM5dxlw\nC707z5dD0Jv7YtvGICVFkGJXMBJWMR/WBq1CwpYQoS6P8URrMeEM5VXbimLvlECC4MYi01C5UzJj\nQV0gYRXzYxGRNiNAwB8VLo/1pQYCQusqf0WYPGmE1AhzJT8bSGV0EeY20Wqq2MbKjkXQbfelUkga\nAJCwijnjzHsGcqFrx3jyq3Vb3zErxyeKash3V9u8PjqiTp+/zfq6CKjL5y4E2xnYN0MkrGLerAlQ\nIAXQXPcJTSzXmpIiCHW7U/Otsfb62tH0qemrs9zAD8eYBNMQLUW/ICSsYv6EotL6eiwlEIroUgQs\nFAF3ThEk2vraXvfBZzekKDXPc6j+DUbCKjaDlTylI0p1il0PUU3qqgf2BdvRI38a2jdFBOr0xdXN\nTxD/5PLLE2kJq9gsQnlG1/auohqKMkPHufa16foDfgGvb/O1N3Q+2gr1kHQtu0BRBSSsYhPxRaeh\nqNV3bJtBoDbd+lDdqcLtbWss35ogkn2EtE2KoVBh7IuEVWwmi9RA6qBMm6i1675FPYt9seNX2tIx\nt5qj67+s35Fvddl0Lr9cUZawis0mJrArdh1F1Vtvh7ys05+RBHMK4Y2lOHzr0XLLFVVAwiq2BdeF\nHItmk4WvRUQZ6t776gnt65NTHZO2kWoz792n7AKQsIrtoU3uMVW0fN38rl332L56PaGy2uxPtVk7\nptHV9/kYEkLXvlhbZoCEVWwXzdTAyvZAuiBFdEP7uuRP29YDpKUpxsL7IxISyhY+zkRUgUzCSvJi\nkneSPFu9H3LYXEDyT0j+OckHSP5SDl/FhtJWYEP7hxDM0L5UkRxTeHwDUSkDVM7yNiun2iRXxHoC\nwN1mdgWAu6v1Jv8E4DVm9lIALwNwNclXTeij2AYGEdgWg1mubV3SDqFtU1MX11QB3NBIdUEuYT0O\n4NZq+VYAb2ga2AH/UK0+o3pt3v/IiDLoI7DOsiaOVn2kzmLoi2+Qr3e58xNVIJ+wXmpmj1bLjwG4\n1GVE8jyS9+Lg77HvNLNPTuWg2FJ8Artm0yJardu5yujk48gj/64ufmiwyuVDil3Qh3mKKjCisJK8\ni+T9jtfxup0d/JuhMxI1s381s5cBOArglSSvDNR3I8k9knv7+/uDtkVsIaHpWT6bNnYhe6C9EPeO\nbBOPXen2J/6QtJ5t8OSsRRUYUVjN7Cozu9Lxuh3AF0geAYDq/fFIWX8H4KMArg7YnDKzXTPbPXz4\n8JBNEdtKG+EcMh+aS1Ta5Ed956ZvPTMX1AW5UgF3ALihWr4BwO1NA5KHSV5ULT8LwOsA/MVkHgqx\nIFVgfcf2Ed2ck/1DI/59UhvOds4/Sq2TS1hvBvA6kmcBXFWtg+TzSZ6ubI4A+CjJ+wDcg4Mc64ey\neCsEEB/57jrJPlbfUMSmRqXkVUPHdvZrcwR1wfk5KjWzJwC81rH9rwFcWy3fB+DlE7smRDq+XOOY\nU6VyTfzvWm9sFsSGojuvhOjLkAIxpXC2iVCbd6I5ywsc17TbYFEFJKxCDIMvTZB7Qv/Q3f8+t6du\ngaAukLAKMQbL+bAZhCQmjGP4FCtzSwR1gYRViLFpe7tnsKzAg06SR/CbE/09UetyOWFqVfBusu0S\nVUDCKsT09BWb+qDZEHNHQ2U1UxdJNwdkjNYLQcIqRE7GEqBYVJqaW3VFu876JKZ1sky3EkI4GDJl\nkFSfYxpV8rMPJKAheHCr/maxu7tre3t7ud0QYlycwthy3betvk8sIXnGzHajdhJWIYRII1VYlWMV\nQoiBkbAKIcTASFiFEGJgJKxCCDEwElYhhBgYCasQQgyMhFUIIQZmI+exktwH8FctDrkEwN+M5M7U\nqC3lsSntANSWrzaz6J/qbaSwtoXkXsqk3zmgtpTHprQDUFtSUSpACCEGRsIqhBADI2E94FRuBwZE\nbSmPTWkHoLYkoRyrEEIMjCJWIYQYmK0UVpIXk7yT5Nnq/VDA9jySf0byQ1P6mEpKW0heTvKjJD9N\n8gGSb8nhqwuSV5N8iOTDJE849pPkr1X77yP5ihx+ppDQlh+o2vApkh8j+dIcfqYQa0vN7ptJPk3y\n+6b0rw0pbSH5apL3VtfHH/Wu1My27gXgHQBOVMsnAPxywPYnAdwG4EO5/e7aFgBHALyiWv5KAJ8B\ncKwA388D8JcAvgbAMwH8edMvANcC+H0ABPAqAJ/M7XePtnwrgEPV8jVzbkvN7g8AnAbwfbn97vG5\nXATg0wBeUK1/Vd96tzJiBXAcwK3V8q0A3uAyInkUwHUA3jORX12ItsXMHjWzP62W/x7AgwAum8xD\nP68E8LCZfdbM/hnA+3HQnjrHAbzPDvgEgItIHpna0QSibTGzj5nZF6vVTwA4OrGPqaR8LgDwYwB+\nF8DjUzrXkpS2XA/gg2b2/wHAzHq3Z1uF9VIze7RafgzApR67dwL4WQBfnsSrbqS2BQBAcgfAywF8\ncly3krgMwOdq6+ewLvgpNiXQ1s834SASL5FoW0heBuC7Abx7Qr+6kPK5vBjAIZJ/SPIMyR/qW+nG\n/pkgybsAPM+x6631FTMzkmtTI0h+J4DHzewMyVeP42UafdtSK+fZOIgwfsLMvjSslyIVkt+BA2H9\nT7l96cE7AfycmX2ZZG5f+nI+gG8C8FoAzwLwcZKfMLPP9ClwIzGzq3z7SH6B5BEze7TqVrpC/28D\n8F0krwVwAYDnkPwNM/vPI7nsZYC2gOQzcCCqv2lmHxzJ1bZ8HsDltfWj1ba2NiWQ5CfJb8RBauka\nM3tiIt/aktKWXQDvr0T1EgDXknzazH5vGheTSWnLOQBPmNk/AvhHkn8M4KU4GIvoRu7kcqaE9q9g\ndcDnHRH7V6PcwatoW3Aw8PM+AO/M7W/Dr/MBfBbAC/HvAwvf0LC5DquDV3+S2+8ebXkBgIcBfGtu\nf/u2pWH/6yh38Crlc/l6AHdXtv8BwP0AruxVb+6GZzrZz61O5FkAdwG4uNr+fACnHfYlC2u0LTjo\nchqA+wDcW72uze175du1OIgM/hLAW6ttNwG4qVomgHdV+z8FYDe3zz3a8h4AX6x9Bnu5fe7aloZt\nscKa2hYAP4ODmQH34yBV1qtO3XklhBADs62zAoQQYjQkrEIIMTASViGEGBgJqxBCDIyEVQghBkbC\nKoQQAyNhFUKIgZGwiq2jeobofSQvIPkfq2dwXpnbL7E56AYBsZWQ/O84eAbEswCcM7O3Z3ZJbBAS\nVrGVkHwmgHsAPIWDe/f/NbNLYoNQKkBsK88F8Gwc/KPCBZl9ERuGIlaxlZC8AwdPk38hgCNm9ubM\nLokNYmOXyIjkAAAAV0lEQVSfxyqEj+oJ8f9iZreRPA/Ax0i+xsz+ILdvYjNQxCqEEAOjHKsQQgyM\nhFUIIQZGwiqEEAMjYRVCiIGRsAohxMBIWIUQYmAkrEIIMTASViGEGJh/A6OTjnSTFsw/AAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fdc79210250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = [0.6, 0.4]\n",
    "x0 = [[-0.4, 0.0, 0.0], [0.6, 0.0, 0.0]]\n",
    "v0 = [[0., -0.2, 0.], [0., 0.3, 0.]]\n",
    "\n",
    "xs = integrate(m, x0, v0, 100, 0.01)\n",
    "plot_xy(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**dt = 0.001, tmax = 100**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVYAAACvCAYAAACmRfbmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEWpJREFUeJzt3X/sXfVdx/HnyzLClGF/CqXt1y9bOhVxi/AVl0ki0xlp\nWSzqYhjTYbKk6ZQ5s6g0sjj+cFrRGLYER7qyrMQhq0pcA13IYDKMjNlv58avBuhqKyUFOqiFLEGo\nvP3jni87vdz7vefee+75nHPu65E03B8f7vl87vec132fzz33HEUEZmZWnh9K3QEzs7ZxsJqZlczB\namZWMgermVnJHKxmZiVzsJqZlczBamZWMgermVnJHKxmZiU7LXUHJmHlypUxOzubuhtm1jL79u37\nXkSsGtSulcE6OzvL/Px86m6YWctIOlyknacCzMxK5mA1MyuZg9XMrGStnGM1q53rfzTBMk9Uv0wD\nHKxm5SkanpMKvO7lL9Yfh+5EOVjNhjUoQFOF1qDl5vvdbwwO3FI4WM0GaUsI9evvYoHbtDHWRNJg\nlXQZ8GlgCbAjIrZ1Pf9B4FpAwEvARyLiO5V31KZLW4K0qF7jWngPHLQjSRaskpYANwG/AhwB9kra\nHRGP5Zr9F/CLEXFc0gZgO/Dz1ffWWq9XmE5ziOTH3q+ineb3Z4CUFevFwIGIOAgg6XZgE/B6sEbE\nA7n2DwJrK+2htZursWK635de1azfu1OkDNY1wFO5+0dYvBr9MPCVifbI2s9hOr5e1axD9hSN+PJK\n0nvoBOsli7TZDGwGmJmZqahn1hgO1MlYLGSn+D1OGaxPA+ty99dmj51C0juAHcCGiHi+34tFxHY6\nc7DMzc1FuV21xnIlVZ2F99cBmzRY9wLrJZ1HJ1CvBK7KN5A0A9wB/E5EPFF9F62RXJ2m1V3FTuGH\nW7JgjYiTkq4B7qZzuNXnI+JRSVuy528G/gxYAfydJICTETGXqs9Wcw7U+pnSKlYR7dtrnpubC5+P\ndYpMYUXUWA3/8JO0r0hx57NbWbM5VJvl+hP9j5FtkUYcFWD2Bg7UZstPEbRwesDBas3iQG2X60+0\ncv7VUwHWHPkNrwUbn2VaOD3gitXqz1XqdGjR9IArVqs3V6nTpwXVqyvWhpvdetfEXvvQtssn9toD\nuUqdbg2vXn0ca0MMCtAyQ7DKZfXUwA3JJqhGH7JFj2N1sNZUv3BLWkVSQb8cqtZLTcLVwdqwYO0V\nWKlDtIhS++1QtUESryMO1gYEa1PDdDHdYyo8HoeqFZVwXXGw1jhYRw6fhik0zprs4lnDJApXB2vN\ngnVawrSXvmN3qNo4EoSrg7UmwTrNgdrLwvtx6Izs1LsOVBtHxeHqYE0crA7UAbINYvbl2/ze2Hgq\nDFcHa6JgdaAWkNsQ8u+X3ysbWUXh6mBNEKwOiQL6bACvTxH4fbNRVRCujTjRtaTLJD0u6YCkrT2e\n/0lJ35D0v5L+KEUfi5jdetcpweBw6GORFX/hPcu/l2ZD6b4MTELJglXSEuAmYANwPvABSed3NXsB\n+APgbyruXmGuUgsqUE3kP5QcrjaSmoRryor1YuBARByMiFeA24FN+QYR8VxE7AVeTdHBQVylFjTk\nLpqrV2u6lMG6Bngqd/9I9thIJG2WNC9p/tixY2N3bjHdu/5WwJDzXq5ebWQ1qFpbcz7WiNgeEXMR\nMbdq1aqJLce7/kMac+V2uNpIEodrymB9GliXu782e6y2HKpDKulbWoerjSThj09SButeYL2k8ySd\nDlwJ7E7Yn0U5VEdU0srtcLWRJahakwVrRJwErgHuBvYDuyLiUUlbJG0BkHSOpCPAx4FPSDoi6ayq\n++ovqUYwgZXZ4WpDSzQlkHSONSL2RMTbI+JtEfGp7LGbI+Lm7PYzEbE2Is6KiKXZ7Rer7KM34jFM\nYFfM4WpN0JovrybNleoQJlwdOFxtKAmqVgfrInxI1Rgm/MWB/yZWZw7WPhyqI6p4LstVqxVScdXq\nYF2EQ3VEFR3m4ikBqysHaw/eUEeU4LAWf/hZYRUe1+pg7eIpgGbyh6EVVkEB4GC1ciX4tYs/BK2w\nitZPB2uOq9Vmc9VqdeFgtXIkPv+lPwxtKBNeXx2sXbyBjqEGV1x11WoDVbCeOlgz3iCbzx+KVhcO\n1hxvmGZWBgerja8GF2/L896HpeZgtXLUYH4VvNdhQ5hgQeBgxRWO2dSZcCHgYDUzK1nSYJV0maTH\nJR2QtLXH85L0mez5hyRdOKm+eBfSzMqSLFglLQFuAjYA5wMfkHR+V7MNwPrs32bgs5V20sxsBAOD\nVdJHJS2bwLIvBg5ExMGIeAW4HdjU1WYTcGt0PAgslbR6An0xMytNkYr1bGCvpF3ZrrtKWvYa4Knc\n/SPZY8O2AUDSZknzkuaPHTtWUhfNzIY3MFgj4hN0dsVvAX4XeFLSX0h624T7NpSI2B4RcxExt2rV\nqtTdMbMpVmiONSICeCb7dxJYBvyTpBvGWPbTwLrc/bXZY8O2KYUPuTKzshSZY/2YpH3ADcC/Az8T\nER8BLgJ+c4xl7wXWSzpP0unAlcDurja7gQ9lRwe8CzgREUfHWKaZ2cQVqViXA78REb8aEf8YEa8C\nRMRrwPtGXXBEnASuAe4G9gO7IuJRSVskbcma7QEOAgeAzwG/N+ryFuNDrcysTKcNahARn1zkuf3j\nLDwi9tAJz/xjN+duB/D74yzDpoundKwQn4/VGqFGJ2LxHogVMsGftTpYc1ztjKgmJ2AxqwsHa8ZV\nTvP5g9HqwsFqreIPSKsDB2sXVz1jqNE8q1lfFaynDtYcVztjSDzP6g9EG4rPx1o9b6TN5A9GqwsH\na5eFjdPhOqIE0wH+W1lhFa2fDlYrT4LpgIVQdbVqhVWwnjpYe3DVOqaKq1aHqhVS4XrpYO3D4Vp/\n/ttYXTlYrVwVTQd4CsCGslCtVrR+OlgX4ap1DBPc7XKoWt05WAdwuI5goSqYQLg6VG1oFVer4GAt\nxOE6ggmsxA5VawoHa0EO1xGVVLU6VG0kCapVSBSskpZL+qqkJ7P/9ry8tqTPS3pO0iNV97EXh+uQ\nSlqZHao2kkShCukq1q3AvRGxHrg3u9/LF4DLqupUEflwdcAWNEbV6lC1JkoVrJuAndntncAVvRpF\nxP3AC1V1qqj8Ru5wHWCML7IcqjayhNUqpAvWs3NXW30GODtRP0Z2aNvlnhooashwze8NOFRtaIlD\nFSYYrJLukfRIj3+b8u2yCwZGCcvbLGle0vyxY8fGfbnCPDVQUMGVPP8eOlStqdTJtYoXKj0OXBoR\nRyWtBu6LiJ/o03YWuDMiLij6+nNzczE/P19KX4tyIBSwSCXh989KMeFqVdK+iJgb1C7VVMBu4Ors\n9tXAlxP1ozTdUwOuXnvoMyWQ3+13qNrIajAFsCBVxboC2AXMAIeB34qIFySdC+yIiI1Zu38ALgVW\nAs8Cn4yIWwa9foqKtZsrsEVkG8Dsy7e9/pDfIxtLRaFatGJNEqyTVodgXeCAfaPZrXdx6IyrOndq\nUF1Yw1VYqTpYaxKsC6Y9YLunRg5tu7xWu27WUBWvQw7WmgXrgp4B02IDx+twtVElWHccrDUN1gW9\nvtxqS8gO/eHhcLVhJVpnHKw1D9a8NoTs2JW4w9WKSriuOFgbFKx5/Q7TqlvQTuTDwOFqi8kfppdo\nHXGwNjRYuw06HnbSgbvY8iey7BpsPFZDNVkvHKwtCdZeUvz4oPKK2dWrLahJqIKDtdXBWlTRAK7b\nNMPrHK7TrUaBusDB6mBthxpuXFaBmv7d636uALNirj8x0YsTWg3VNFSHcVrqDpgVcv2Jzgbn6YH2\nakGgLnCwWnPkK1cHbHt074m04G/qqQBrnvyG5+mBZuuuUlsQquCK1ZrK1WuztWi3vxcHqzWbA7ZZ\nWh6oCxys1g4O2HqbkkBd4GC1dukVsPnHrTot/FKqqCTBKmk58CVgFjhE59Isx7varANupXNp7AC2\nR8Snq+2pNVb3sa+uYqszxYG6INU1r24AXoiIbZK2Assi4tquNquB1RHxLUlvAfYBV0TEY4Ne37+8\nsp68wU9Or6MzWvj+Fv3lVaqpgE10LhIIsBO4DzglWCPiKHA0u/2SpP3AGmBgsJr15Cq2XFMSpqNI\nVbH+T0QszW4LOL5wv0/7WeB+4IKIeLFPm83AZoCZmZmLDh8+XHa3rY1cxQ5nysM0+UlYJN0DnNPj\nqeuAnfkglXQ8Ipb1eZ0zga8Dn4qIO4os21MBNpIpD42+/L68LvlUQES8t99zkp6VtDoijmZzqc/1\nafcm4J+BLxYNVbOR9fpF1zRWtP1+zTYNYy9JqjnW3cDVwLbsv1/ubpBNEdwC7I+Iv622ezb1ukOk\nX9D2atski/0kuMnjSizVHOsKYBcwAxymc7jVC5LOBXZExEZJlwD/BjwMvJb9r38aEXsGvb6nAqwS\ng85TUKdgKnJOhTr1t6aSz7Gm5GC1ZEY5Kcw4gVb18qZc8jlWs6lUJLTeMG9b8hm6HJzJOVjNqubg\naz2fj9XMrGQOVjOzkjlYzcxK5mA1MyuZg9XMrGStPI5V0jE6PzwoaiXwvQl1p2oeS/20ZRzgsfx4\nRKwa1KiVwTosSfNFDvptAo+lftoyDvBYivJUgJlZyRysZmYlc7B2bE/dgRJ5LPXTlnGAx1KI51jN\nzErmitXMrGQOVjOzkk1lsEpaLumrkp7M/tvzeltZ2yWS/lPSnVX2sagiY5G0TtK/SnpM0qOSPpai\nr71IukzS45IOZJdC735ekj6TPf+QpAtT9LOIAmP5YDaGhyU9IOmdKfpZxKCx5Nr9nKSTkt5fZf+G\nUWQski6V9O1s+/j62AuNiKn7B9wAbM1ubwX+apG2HwduA+5M3e9RxwKsBi7Mbr8FeAI4vwZ9XwJ8\nF3grcDrwne5+ARuBrwAC3gV8M3W/xxjLu4Fl2e0NTR5Lrt3XgD3A+1P3e4y/y1LgMWAmu/9j4y53\nKitWYBOwM7u9E7iiVyNJa4HLgR0V9WsUA8cSEUcj4lvZ7ZeA/cCaynrY38XAgYg4GBGvALfTGU/e\nJuDW6HgQWJpdgLJuBo4lIh6IiOPZ3QeBtRX3sagifxeAj9K52GfPi4HWRJGxXAXcERH/DRARY49n\nWoP17Ig4mt1+Bji7T7sbgT/hB9fcqqOiYwFA0izws8A3J9utQtYAT+XuH+GNgV+kTR0M288P06nE\n62jgWCStAX4d+GyF/RpFkb/L24Flku6TtE/Sh8ZdaGuvICDpHuCcHk9dl78TESHpDcecSXof8FxE\n7JN06WR6Wcy4Y8m9zpl0Kow/jIgXy+2lFSXpPXSC9ZLUfRnDjcC1EfFa54LKjXYacBHwy8CbgW9I\nejAinhjnBVspIt7b7zlJz0paHRFHs93KXqX/LwC/JmkjcAZwlqS/j4jfnlCX+yphLEh6E51Q/WJE\n3DGhrg7raWBd7v7a7LFh29RBoX5KegedqaUNEfF8RX0bVpGxzAG3Z6G6Etgo6WRE/Es1XSysyFiO\nAM9HxPeB70u6H3gnne8iRpN6cjnRhPZfc+oXPjcMaH8p9f3yauBY6HzxcytwY+r+dvXrNOAgcB4/\n+GLhp7vaXM6pX179R+p+jzGWGeAA8O7U/R13LF3tv0B9v7wq8nf5KeDerO0PA48AF4y13NQDT/Rm\nr8jeyCeBe4Dl2ePnAnt6tK9zsA4cC51dzgAeAr6d/duYuu9Z3zbSqQy+C1yXPbYF2JLdFnBT9vzD\nwFzqPo8xlh3A8dzfYD51n0cdS1fb2gZr0bEAf0znyIBH6EyVjbVM/6TVzKxk03pUgJnZxDhYzcxK\n5mA1MyuZg9XMrGQOVjOzkjlYzcxK5mA1MyuZg9WmTnYO0YcknSHpR7JzcF6Qul/WHv6BgE0lSX9O\n5xwQbwaORMRfJu6StYiD1aaSpNOBvcDLdH67/3+Ju2Qt4qkAm1YrgDPpXFHhjMR9sZZxxWpTSdJu\nOmeTPw9YHRHXJO6StUhrz8dq1k92hvhXI+I2SUuAByT9UkR8LXXfrB1csZqZlcxzrGZmJXOwmpmV\nzMFqZlYyB6uZWckcrGZmJXOwmpmVzMFqZlay/wcIGpIOTeoY4wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fdc82160510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = [0.6, 0.4]\n",
    "x0 = [[-0.4, 0.0, 0.0], [0.6, 0.0, 0.0]]\n",
    "v0 = [[0., -0.2, 0.], [0., 0.3, 0.]]\n",
    "\n",
    "xs = integrate(m, x0, v0, 100, 0.001)\n",
    "plot_xy(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Test for three-body problem\n",
    "\n",
    "**initial condition:**\n",
    "\n",
    "$$m_1 = 3, m_2 = 4,\\text{ and } m_3 = 5$$\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray*}\n",
    "\\textbf{r}_1 &=& (1.0, 3.0, 0.0)\\\\\n",
    "\\textbf{r}_2 &=& (-2.0, -1.0, 0.0)\\\\\n",
    "\\textbf{r}_3 &=& (1.0, -1.0, 0.0)\\\\\n",
    "\\textbf{v}_{1,2,3} &=& (0, 0, 0)\n",
    "\\end{eqnarray*}\n",
    "$$\n",
    "\n",
    "use G = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**dt = 0.001, tmax = 5.0**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALYAAAFACAYAAAARTIa0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADwlJREFUeJzt3X+o3fddx/Hna62zspWRmhBj23hTyYRswYjX4o9Zpvth\n1wpZB45W2CIMUqEtE0SJP2D5w2koq/tLqq0tyx92pazGhk1alioOmWhuRmiShpqs3tKEtEkddBNc\nt7Rv/7jfG09vzr333Jzv9/v58X094JJ7zznJfROefHjfc+89X0UEZrV5R+oBzLrgsK1KDtuq5LCt\nSg7bquSwrUoO26rksK1KDtuqdHXqAdqwfv36mJmZST2G9eDIkSOvRcSG1R5XRdgzMzPMzc2lHsN6\nIOmlSR7nVcSq5LCtSg7bquSwrUoO26rksK1KDtuq5LCtSg7bquSwLU973zPVX08WtqQbJf2zpOcl\nnZD02eb26yR9XdKp5s91qWa0cqU8sS8Cvx8R24BfBO6RtA3YAzwbEVuBZ5uPbUimPK0hYdgRcS4i\nvtW8/z3gJHA9sBPY3zxsP/DxNBNaUntfn+qvZ7FjS5oBfg74d2BjRJxr7noF2JhoLCtY8rAlvRt4\nEvi9iPju6H2x8DJVY1+qStJuSXOS5i5cuNDDpNaLFtYQSBy2pB9hIeq/i4i/b25+VdKm5v5NwPlx\nfzciHoqI2YiY3bBh1Z87t5JMuYZA2mdFBDwCnIyIvxy56yCwq3l/F/BU37NZ+VL+Bs2vAJ8Cjkk6\n2tz2x8A+4AlJnwFeAj6ZaD4rWLKwI+JfAS1z94f6nMUy0dJ+DRl88Wj2Ni3s1+CwrVIO2/LQ4hoC\nDtty0tIaAg7bKuWwLb2W1xBw2JaLFtcQcNhWKYdtaXWwhoDDthy0vIaAw7ZKOWyrksO2dPa+p5M1\nBBy2VcphW5UctqXR0dN8ixy2pdPRfg0O2yrlsK1/Ha8h4LAtlQ7XEHDYVimHbf3qYQ0Bh20pdLyG\ngMO2Sjls609Pawg4bOtbD2sIOGyrVOqXEX5U0nlJx0du2yvprKSjzdttKWe0MqU+sb8E3Drm9i9G\nxI7m7R97nsm60ON+DYnDjohvAN9JOYP1qKf9GtKf2Mu5T9Jzzaoy9nJ4vlSHrSTHsB8EbgJ2AOeA\nB8Y9yJfqKEjPawhkGHZEvBoRb0bEW8DDwM2pZ7IW9LiGQIZhL15YqXEHcHy5x5otJ+U1aJD0ZeCD\nwHpJZ4DPAR+UtIOFy+DNA3cnG9Cml2ANgcRhR8RdY25+pPdBrFs9ryGQ4Spi1gaHbd1JtIaAw7au\nJVhDwGFbpRy2dSPhGgIO27qUaA0Bh22VcthWJYdt7Uu8X4PDtq4k3K/BYVulHLa1K4M1BBy2dSHx\nGgIO2yrlsK09mawh4LCtbRmsIeCwrVIO29qR0RoCDtvalMkaAg7bKuWwrUoO26aX2X4NDtvaktF+\nDQ7bKuWwbToZriHgsK0Nma0hkOelOq6T9HVJp5o/x74+ttlKUp/YX+LyS3XsAZ6NiK3As83HlqNM\n1xDI81IdO4H9zfv7gY/3OpStTYZrCKQ/scfZGBHnmvdfATamHMbKlGPYl0REsPA62ZfxNWgSy3gN\ngTzDfnXxqgbNn+fHPcjXoMlApmsI5Bn2QWBX8/4u4KmEs1ihUj/d92Xg34CfkXRG0meAfcBHJJ0C\nPtx8bDnJfA2BPC/VAfChXgextct4DYE8VxGzqTlsq5LDtrUpYL8Gh21XIvP9Ghy2Vcph2+QKWUPA\nYdtaFbCGgMO2Sjlsm0xBawg4bFuLQtYQcNhWKYdtqytsDQGHbZMqaA0Bh22VcthWJYdtKytwvwaH\nbZMobL8Gh22Vcti2vELXEHDYtpoC1xBw2FYph23jFbyGgMO2lRS6hoDDtko5bLtc4WsIOGxbTsFr\nCDhsq1TS1+5biaR54HvAm8DFiJhNO9FAVLCGQMZhN34tIl5LPcTgFL6GgFcRq1TOYQdwSNIRSbuX\n3ulLddhKcg77AxGxA/gYcI+kW0bv9KU6OlDJfg0Zhx0RZ5s/zwMHgJvTTjQQFezXkGnYkt4l6drF\n94GPAsdX/ltm/y/XZ0U2AgckwcKMj0XE02lHqlxFawhkGnZEvAj8bOo5BqeSNQQyXUXMpuWwrbo1\nBBy2LapoDQGHbZVy2ENX4RoCDtugujUEHLZVatWwJd0naV0fw1j3ZvZ8LfUIvZjkxN4IHJb0hKRb\n1Xw70CpQ6X4NE4QdEX8KbAUeAX4HOCXpzyX9dMezWR8q3K9hwh07IgJ4pXm7CKwDviLp/g5ns5YN\nZQ2BCX5WRNJngU8DrwF/C/xBRPxQ0juAU8AfdjuitWl+3+0L71S8hsBkPwR1HfCJiHhp9MaIeEvS\nb3YzlvWi0jUEJgg7Ij63wn0n2x3HrB1+Hnsg3rZfV76GgMMelEv7NVS9hoDDtko57KEZwBoCDntQ\nLu3Zla8h4LAHY37f7cxf89upx+iNwx6Yme8/lnqEXjjsoRjZrYfwrXWHPSR7X7/0lF/tcTvsIVjy\nTMgQ4nbYQ7HkmZDa43bYtVvheeua48427Oa3dV6QdFrSntTzFGkx6hWet6417izDlnQV8FcsvDb2\nNuAuSdvSTlWoCb4ZU2PcWYbNwmthn46IFyPiB8DjwM7EM5Vljd86ry3uXMO+Hnh55OMzzW22Fmv8\n1nlNceca9qp8DZoVTPGDTrXEnWvYZ4EbRz6+obntEl+DZhVT/KBTDXHnGvZhYKukLZLeCdwJHEw8\nUxla+rHU0uPOMuyIuAjcCzwDnASeiIgTaafK3/b929m+ZTPbt2xu5d8rOW4tvGRI2WZnZ2Nubi71\nGP1bcjpv37KZY7uOsX3/9sseemzXsSv+NIthv+1XyxKRdGSSy4877Nyttlo0u/RizOMCXhr6lUSe\nS9wOuwYTfOcQVo563GOv9PTOIe5Jw85yx7bG3tcnfnZj0linWUlK2rkdduHG7dNdKiVuh12wtawg\nbSohbodduL6jXpR73A67UH2vIOPkHLfDLliq03pUrnE77ALlcFqPyjFuh12YVF8wria3uB12gXKL\nelFOcTvsguS2goyT+lvuixx2IXJdQcaZ33d78lPbYRekhKhHpYzbYReghBVkqdT7tsMuRGmnNaSN\n22FnrsTTelSquB12xkr6gnElKeJ22JkrPepFfT8N6LAzVfoKMk6fcTvsDNWygqTksDPlqKfjsDNT\n4wqSgsPOkE/r6TnsjPi0bo/DzoS/YGyXw86Io25PdmFL2ivprKSjzdttqWfqmleQ9l2deoBlfDEi\nvpB6iD54BelGdif2EDnq9uUa9n2SnpP0qKR14x5Qw6U6vIJ0J0nYkg5JOj7mbSfwIHATsAM4Bzww\n7t+o5VIdPq27kWTHjogPT/I4SQ8DX+14nCR8Wncru1VE0qaRD+8AjqeapSv+grF7OT4rcr+kHUAA\n88DdacfphqPuVnZhR8SnUs/QJa8g/chuFRkCn9bdc9g9mub6L7Y2DrsnXkH65bB75NO6Pw67Bz6t\n++ewO+bnrNNw2D1w1P1z2B3yCpKOw+6YT+s0HHZHfFqn5bA74C8Y03PYHXHUaTnslnkFyYPDbpFX\nkHw47JY56jw47JZ4BcmLw26RT+t8OOwW+LTOj8Oekr9gzJPDboGjzo/DnoJXkHw57CvkFSRvDnsK\njjpfDvsKeAXJn8O+Qj6t8+aw18indRlSvYzwb0k6IektSbNL7vsjSaclvSDpN1LMtxx/wViOVK/d\ndxz4BPA3ozdK2gbcCbwP+EngkKT3RsSb/Y84nqMuQ5ITOyJORsQLY+7aCTweEW9ExH8Bp4Gb+51u\nPK8gZcltx74eeHnk4zPNbUl5BSlPZ6uIpEPAT4y5608i4qkW/v3dwG6AzZs3T/vPrcpRl6WzsCe9\nHMcSZ4EbRz6+oblt3L//EPAQwOzsbFzB55qIV5Ay5baKHATulPSjkrYAW4H/SDyTT+sCpXq67w5J\nZ4BfAr4m6RmAiDgBPAE8DzwN3JPyGRGf1uVKddWwA8CBZe77PPD5fie6nL9gLFtuq0hWHHW5HPYY\nXkHK57CX4dO6bA57CV8AqQ4Oe4RXkHo47CV8WtfBYTd8WtfFYePnrGvksBuOui6DD9srSJ0GHzb4\ntK7RoMP2aV2vwYbtLxjrNtiwwVHXbJBhewWp3+DC9goyDIMLGxz1EAwqbK8gwzGosMGn9VAMJmyf\n1sMyiLD9BePwDCJscNRDU33YXkGGqeqwvYIMV9Vhg6MeqqrDdtTDVXXYNlxZXYNG0oyk/5V0tHn7\n6xTzWfmyugZN49sRsaPneawyqV5t9SSApBSf3gYgxx17S7OG/IukX13uQZJ2S5qTNHfhwoU+57MC\n5HYNmnPA5oj4b0k/D/yDpPdFxHeXPrCvS3VYmbK6Bk1EvAG80bx/RNK3gfcCcy2PZ5XLahWRtEHS\nVc37N7FwDZoX005lJcrqGjTALcBzko4CXwF+NyK+k2JGK1tW16CJiCeBJ/ufyGqT1Spi1haHbVVS\nRPnPlEm6ALyUeIz1wGuJZ1iU0yzQ7jw/FREbVntQFWHnQNJcRMyu/sju5TQLpJnHq4hVyWFblRx2\nex5KPcCInGaBBPN4x7Yq+cS2Kjlsq5LDbomkvZLOjvxa220JZrhV0guSTkva0/fnHzPPvKRjzf9H\nrz+h6R27JZL2Av8TEV9I9PmvAv4T+AhwBjgM3BURz6eYp5lpHpiNiN6/WeQTux43A6cj4sWI+AHw\nOLAz8UzJOOx23SfpOUmPSlrX8+e+Hnh55OMzzW0pBXBI0hFJu/v8xA57DSQdknR8zNtO4EHgJmAH\nC7/i9kDSYfPwgeYVBz4G3CPplr4+caqXXyjSpL/uJulh4Ksdj7PUWeDGkY9vaG5LJiLONn+el3SA\nhXXpG318bp/YLZG0aeTDO1h47ZQ+HQa2Stoi6Z3AncDBnme4RNK7JF27+D7wUXr8P/GJ3Z77Je1g\nYa+cB+7u85NHxEVJ9wLPAFcBj0bEiT5nWGIjcKB57Zirgcci4um+Prmf7rMqeRWxKjlsq5LDtio5\nbKuSw7YqOWyrksO2Kjnswkj6heYHra5pvrt3QtL7U8+VG3+DpkCS/gy4Bvgx4ExE/EXikbLjsAvU\n/CzIYeD7wC9HxJuJR8qOV5Ey/TjwbuBaFk5uW8IndoEkHWThN2S2AJsi4t7EI2XHP91XGEmfBn4Y\nEY81v+f4TUm/HhH/lHq2nPjEtip5x7YqOWyrksO2Kjlsq5LDtio5bKuSw7Yq/R+LRpMqZxlUMgAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fdc791e9d50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = [3.0, 4.0, 5.0]\n",
    "x0 = [[1.0, 3.0, 0.0], [-2.0, -1.0, 0.0], [1.0, -1.0, 0.0]]\n",
    "v0 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]\n",
    "\n",
    "xs = integrate(m, x0, v0, 5.0, 0.001)\n",
    "plot_xy(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**dt = 0.0001, tmax = 5.0**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAEzCAYAAABaGjpLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEKRJREFUeJzt3X+IZWd9x/HPJ5uUSCPadocqSbazf0gxZEpDp2mt/5RU\ny9ZIRUHQgsa2sP80JYKtbEmLKRS6UBAKCu1igoGmSqgWQ1fZrhgJ0hozG1I3m9U2tdm6omSDrUZK\nW9d8+8fcu717vTPz3LnnnOfHeb9gYWdndua5JPvm+zzn3nMdEQIA7O2a3AsAgFoQTABIRDABIBHB\nBIBEBBMAEhFMAEhEMAEgEcEEgEQEEwASXZt7Acs4ePBgrK+v514GgMacOXPmhYhY2+vrqgrm+vq6\ntra2ci8DQGNsX0j5OrbkAJCIYAJAIoIJAIkIJgAkIpgAkIhgAkCibMG0fb3tL9n+J9vnbP9xrrUA\nQIqcz8P8H0l3RMT3bF8n6Qu2PxMRX8y4JgDYUbZgxvabCX1v8uF1k1+8wRCAYmU9w7R9wPZTkp6X\ndDoiHs+5HqBE68dO5l4CJrIGMyJ+EBE/K+kmSbfbvnX+a2wftb1le+vSpUvDLxIAJoq4Sh4R/ynp\nUUlHFnzuRERsRsTm2tqer40HgN7kvEq+ZvuVk9+/TNIbJX0l13qAErEdL0vOq+SvlvSg7QPaDvfD\nEfF3GdcDFOm543fmXgImcl4l/7Kk23L9fABYVhFnmABQA4IJFIrzy/IQTKBgnF+WhWACQCKCCQCJ\nCCZQIM4vy0QwgUJxflkeggkAiQgmACQimACQiGACheGCT7kIJlAgLviUiWACQCKCCQCJCCZQEM4v\ny0YwgcJwflkuggkAiQgmACQimEAhOL8sH8EECsL5ZdkIJgAkIpgAkIhgAgXg/LIOBBMoBOeX5SOY\nAJCIYAJAIoIJAIkIJpAZF3zqQTCBAnDBpw4EEwASEUwASEQwgYw4v6wLwQQy4/yyHgQTABIRTCAT\ntuP1IZhARmzH60IwASARwQSARAQTyIDzyzoRTCATzi/rQzABIBHBBIBEBBMYGOeX9SKYQAacX9aJ\nYAJAIoIJAImyBdP2zbYftf2M7XO278m1FmAonF/W7dqMP/uypPdFxJO2Xy7pjO3TEfFMxjUBveP8\nsl7ZJsyI+GZEPDn5/YuSzku6Mdd6AGAvRZxh2l6XdJukxxd87qjtLdtbly5dGnppQGfYjtcvezBt\n3yDpE5LeGxHfnf98RJyIiM2I2FxbWxt+gUCH2I7XLWswbV+n7Vg+FBGfzLkWANhLzqvklnS/pPMR\n8cFc6wCAVDknzNdLepekO2w/Nfn1pozrAXrD+WUbsj2tKCK+IMm5fj4wNM4v65f9og8A1IJgAj1j\nO94OggkMgO14GwgmACQimACQiGACPeL8si0EE+gZ55ftIJgAkIhgAj1hO94eggn0iO14WwgmACQi\nmACQiGACPeD8sk0EE+gJ55ftIZgAkIhgAh1jO94uggn0gO14mwgmACQimECH2I63jWACHWM73i6C\nCQCJCCYAJCKYQEc4v2wfwQQ6xPll2wgmACQimEAH2I6PA8EEOsJ2vH0EEwASEUwASEQwgRVxfjke\nBBPoAOeX40AwASARwQRWwHZ8XAgmsCK24+NBMAEgEcEE9ont+PgQTGAFbMfHhWACQCKCCQCJCCaw\nD5xfjhPBBPaJ88vxIZgAkIhgAktiOz5eBBPYB7bj40QwASBR1mDafsD287afzrkOAEiRe8L8qKQj\nmdcAJOP8ctyyBjMiHpP07ZxrAJbF+eV45Z4w92T7qO0t21uXLl3KvRwAI1Z8MCPiRERsRsTm2tpa\n7uUAGLHigwmUgvNLEExgCZxfjtu1OX+47Y9J+mVJB21flPSBiLg/55rQro0HN5b+O2fvOtvDSlCr\nrMGMiHfm/Plo26JAEkCsImswgT7MhrKrQHJ+CYlgohHz02QfkyTnl9gzmLZ/V9JfRcR/DLAeYCl9\nTJPATlImzJ+U9ITtJyU9IOlURES/ywJ2RyiRw57BjIg/tP1Hkn5V0m9K+pDthyXdHxH/2vcCgVk5\nQsn5JaaSzjAjImx/S9K3JF2W9GOS/sb26Yh4f58LBKT8EyXnl5DSzjDvkfRuSS9I+oik34+I79u+\nRtK/SCKY6NU0lmy9kVvKhPnjkt4WERdm/zAiXrL95n6WBeSfKoF5KWeYH9jlc+e7XQ5QVig5v8Qs\nnoeJopS4/eb8ElMEE0UoaaoEdkIwkV2JUyWwCLd3Q1Ylx5LzS8xjwkQWtWzBOb/ELIKJwZU8VQK7\nYUuOQRFL1IxgYjDEErUjmBgEsUQLCCZ6V2MsuUKORQgmelVjLKe4Qo55BBO9qTmWwCI8rQiL3feK\nJb72Oz/0R8QSLSKYY7dbGBeEcOHfn/seG4cPSSKWaA/BHJtFgUwJ447fb+7vTr7/2X/79/1/T6BQ\nBHMM5iO5SiD3MJ0ur/zcHn9WX7hCjp0QzJbNhnLAcF3Zii9zDloYrpBjEYLZokyhXKjSKRNYhGC2\npKRQTtdQ8ZQJzCOYLSgolGfvOquNBze2t+XEEo0hmDUrKJTzNh7ckA4f4qlFaArBrNU0loWFUtKV\n6XLj8KFqbhQMpCCYtSl4qrxi+lzMmUBuPLhxVTyniChq4ojIvYZkm5ubsbW1lXsZ+RQ8VV6x5BoX\nRbREhL1tts9ExOZeX8eEWYsGYymVGaL1Yyeveh7mTtPxVImPAf0gmDVoNJYlWvQqn92CuCimBLRd\nBLN0NYSohjUuYZlX+czHkYC2jWCWrPQQ1XABamC7BZRw1o9glq7EEBHKZNNIEs42EMxSlfoqmdKn\n3kIRzjbwFhUlKylKszcKLmldlTl719mr4om6MGGWqKTpku13L6avuWfarAvBLFXOOA14w+Exm9+m\nE83ysSUvTc7pcn7bPf01MkPfcZ0tej0IZomGjNQ0kvOhHKlcd1onmnVgSz5GbLmLND3XRLkI5lgQ\nyWpwnlmurMG0fUTSn0s6IOkjEXE853qaQySrw5RZtmzBtH1A0oclvVHSRUlP2H4kIp7JtaYmEMkm\nMGWWKeeEebukZyPia5Jk++OS3iKJYC5j0VV1Ilk1psxy5QzmjZK+PvPxRUm/MP9Fto9KOipJhw4d\nGmZlJSOQo8GUWZ7iL/pExAlJJ6TtO65nXs4wZt/Lm0BmMX8T4aExZZYpZzC/IenmmY9vmvzZeF31\nMkRekpjLc8fvHPzJ66hDzmA+Iek1tg9rO5TvkPQbGdczrJ1e0XPfd7jJBVCobMGMiMu275Z0SttP\nK3ogIs7lWk9vdnup405BnI0mgGJkPcOMiE9L+nTONXRmP2FM+Z5MmUAxir/oU4TUaa/LuDFlZpf7\nwg/K034wu4pOrkmPKTMLLvxgkREEs+LYTKdMogkUgdu7lW6352Oid0yZmEUwa0A0s+D8EvMIZi2I\nZjY5pkxe5VMmglkTojm43C+PRFkIZm2IJpANwazRbDQJ5yCG3JazHS8XwazV7JuVEc1e5diWsx0v\nE8GsHdPmYIaYMpkuy0YwW8C02bshp0ymy3IRzJYwbfauzymT6bJ8BLM189Mm4ezMdMrsM5pMl2Uj\nmK0inL3oa2vOdFkHgtk6wtmLLqfMaSyZLsvniHreV2xzczO2trZyL6NuvG95J6bBXHXiJJZlsH0m\nIjb3+jomzLGZTpxMnSvpYmtOLOtDMMdsUTiJ51L2uzUnlnViS46r8fa+S9nP1pxYlid1S04wsRhn\nncmWiSaxLBPBRHeI555Sokksy5UazPbf0wermw3k/Dkn8ZS095umEcs2EEwsh3juav6teWefkE4s\n68eWHN1g2y7p6q05U2U9OMNEPiOPJ1NlfTjDRD4j3bbPhvLF88czrgR9YcLEcBqdPHeaKLt6+ST6\nx5YcZWsgnilb79kr54SzXAQT9ahoyz5/G7bUM0qmzbIRTNRnEs6Nw4dW/lZdXmzZbyTnEc1yEUxU\nqYun4ux2M97U79tVJOcRzTIRTNTnvldcNV12/ZScRSE9e9fZHQPb11OCiGZ5CCbKt+hWcpMzzK4j\nlvoWEEM9b5KLQWUhmCjLTvfZXOIiz6rve7NXDHeaQPvEtFkGgok8drsBceFXwOf1dY45j2jmRzDR\nn73uyl5ZGFP0HU+26HkRTKxmhFFM1Wc8mTbzIJjYG1FcWR/xZNocHsHE/2voXLFkXceTaXM4BHNs\nmBaL0vWrgyTC2SeC2SKiWKUu4kk4+0Uwa0UUm7bqzYXZpveDYJaMKI7eKlMn02b3CGZuRBGJ9jt1\nEs7uFB1M22+XdJ+k10q6PSKSKlhkMLkCjY6seq9NiXDuV+nBfK2klyT9paTfKz6YTIsY2H6mTsK5\nf0UH88oPtz+vUoPZwFsooH77uT8o4VxeM8G0fVTSUUk6dOjQz124cGGg1QF1I5zpsgfT9mclvWrB\np+6NiE9NvubzKnXCBBpBOPeW/X3JI+INfX1vAOmmkVw/dpLnca7omtwLADCM547fuTCeSJclmLbf\navuipNdJOmn7VI51AGNEOPePJ64DIzcfzDFu17Nf9OkDwQT6NdYLRAQTwL6NLZwEE8DKxrJdJ5gA\nOtVyPAkmgN60tmUnmAB618rUSTABDKrmeBJMANnUFk+CCaAIi15JVFpACSaAIvUd0PVjJ5f+ftnv\nVgQAi8zHbKfXsy8TvaFeD8+ECaA4+wngKlMqEyaAapV2xjnF/TABIBHBBIBEBBMAEhFMAEhEMAEg\nEcEEgEQEEwASEUwASEQwASARwQSARAQTABJVdfMN25ckXRjoxx2U9MJAP2tIrT4uqd3HxuPq309F\nxNpeX1RVMIdkeyvl7iW1afVxSe0+Nh5XOdiSA0AiggkAiQjmzk7kXkBPWn1cUruPjcdVCM4wASAR\nEyYAJCKYAJCIYO7C9p/Z/ortL9v+W9uvzL2mLth+u+1ztl+yXdXTOhaxfcT2V20/a/tY7vV0xfYD\ntp+3/XTutXTJ9s22H7X9zOT/w3tyrykVwdzdaUm3RsTPSPpnSX+QeT1deVrS2yQ9lnshq7J9QNKH\nJf2apFskvdP2LXlX1ZmPSjqSexE9uCzpfRFxi6RflPQ7tfw3I5i7iIi/j4jLkw+/KOmmnOvpSkSc\nj4iv5l5HR26X9GxEfC0i/lfSxyW9JfOaOhERj0n6du51dC0ivhkRT05+/6Kk85JuzLuqNAQz3W9J\n+kzuReCH3Cjp6zMfX1Ql//gg2V6XdJukx/OuJM3o35fc9mclvWrBp+6NiE9NvuZebW8jHhpybatI\neVxATrZvkPQJSe+NiO/mXk+K0QczIt6w2+dtv0fSmyX9SlT0pNW9HldDviHp5pmPb5r8GQpm+zpt\nx/KhiPhk7vWkYku+C9tHJL1f0q9HxH/lXg8WekLSa2wftv0jkt4h6ZHMa8IubFvS/ZLOR8QHc69n\nGQRzdx+S9HJJp20/Zfsvci+oC7bfavuipNdJOmn7VO417dfkotzdkk5p++LBwxFxLu+qumH7Y5L+\nUdJP275o+7dzr6kjr5f0Lkl3TP5dPWX7TbkXlYKXRgJAIiZMAEhEMAEgEcEEgEQEEwASEUwASEQw\nASARwQSARAQTTbH985P7l15v+0cn91u8Nfe60AaeuI7m2P4TSddLepmkixHxp5mXhEYQTDRn8pry\nJyT9t6RfiogfZF4SGsGWHC36CUk3aPs+ANdnXgsawoSJ5th+RNt3Xj8s6dURcXfmJaERo78fJtpi\n+92Svh8Rfz15v59/sH1HRHwu99pQPyZMAEjEGSYAJCKYAJCIYAJAIoIJAIkIJgAkIpgAkIhgAkCi\n/wMLQUv9OsGxdwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fdc792cb9d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = [3.0, 4.0, 5.0]\n",
    "x0 = [[1.0, 3.0, 0.0], [-2.0, -1.0, 0.0], [1.0, -1.0, 0.0]]\n",
    "v0 = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]\n",
    "\n",
    "xs = integrate(m, x0, v0, 5.0, 0.0001)\n",
    "plot_xy(xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}