rk2_simple_harmonic_oscillator.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "iterations = 100\n",
      "#global y\n",
      "y = 1.\n",
      "#global dy\n",
      "dy = 0.\n",
      "h = 0.1 # step size\n",
      "w = 1. # omega or the time period\n",
      "#k1, k2, K1, K2 = 0, 0, 0, 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 72
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rk2ylist = []\n",
      "rk2dylist = []\n",
      "rk2ylist.append(y)\n",
      "rk2dylist.append(dy)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 79
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rk2y = 1.\n",
      "rk2dy= 0.\n",
      "for i in xrange(iterations):\n",
      "    #k_update()\n",
      "    rk2k1 = rk2dy # y update\n",
      "    rk2K1 = -rk2y # dy update\n",
      "    rk2k2 = rk2dy + rk2K1*h # y update\n",
      "    rk2K2 = -(rk2y+h)# dy update    \n",
      "    #print k1, k2, K1, K2\n",
      "    rk2y += h*(rk2k1+rk2k2)/2.\n",
      "    rk2dy += h*(rk2K1+rk2K2)/2.\n",
      "    #system_update()\n",
      "    #y += h*(k1+k2)/2.\n",
      "    #dy += h*(K1+K2)/2.\n",
      "    rk2ylist.append(rk2y)\n",
      "    rk2dylist.append(rk2dy)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 80
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(ylist,'o',dylist,'x')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 71,
       "text": [
        "[<matplotlib.lines.Line2D at 0xb153f22c>,\n",
        " <matplotlib.lines.Line2D at 0xb153f66c>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX90FPW5/98roSS2QCJXfphEgsuPSAIxiEXlQGJpspoA\n19oeQYtwMe1VoYH2CPi9JhziFSICvTVBW2xPesu3vRpsa4VvFsJiuZt4QgJX8KIIR2zcQBIIKu6W\n1oZA8Pn+scxmdnZ2k+zsj5n5PK9zOOzsTGY+Ozv7PJ/P87w/z8dCRASGYRhGOG6IdwMYhmGY+MAO\ngGEYRlDYATAMwwgKOwCGYRhBYQfAMAwjKOwAGIZhBEWzA3j88ccxZswYTJs2TXW/0+nEyJEjkZub\ni9zcXGzcuFHrJRmGYZgIkKD1BMuXL0dpaSmWLl0a9Ji8vDzs2bNH66UYhmGYCKJ5BDBnzhykpKSE\nPIbnmjEMw+iPqOcALBYLDh06hJycHBQVFeHkyZPRviTDMAwzADSHgPpjxowZaG9vx4033oh9+/bh\nwQcfxOnTp6N9WYZhGKYfou4Ahg8f7nv9wAMPYMWKFfjiiy9w0003+R1nsVii3RSGYRhTEm6YPeoh\noAsXLvgad+TIERBRgPGXICL+R4QNGzbEvQ16+cf3gu8F34vQ/7SgeQTwyCOPoKGhAZ9//jnS09Px\n3HPP4erVqwCAJ554An/4wx/wi1/8AgkJCbjxxhtRW1ur9ZIMwzBMBNDsAF5//fWQ+1euXImVK1dq\nvQzDMAwTYXgmsA7Jz8+PdxN0A9+LPvhe9MH3IjJYSGsQKUJYLBbN8SyGYRjR0GI7eQTAMAwjKOwA\nGIZhBIUdAMMwjKCwA2AYhhEUdgAMwzCCwg6AYRhGUNgBMAzDCAo7AIZhGEFhB8AwDCMo7AAYhmEE\nhR0AwzCMoLADYBiGERR2AAzDMILCDoBhGEZQ2AEwDMMICjsAhmEYQWEHwDAMIyjsABiGYQSFHQDD\nMIygsANgGIYRFHYADMMwgsIOgGEYRlDYATAMwwhKQrwbwDCMGNjtjaiudqCnJwGXLnUA+BpGjBiN\nYcN6cc89t6C5+ZzqvlWrClFcPDfezTclFiKieDcCACwWC3TSFIZhIozd3ojVq/ejtXUTgEYA+wFs\nur63EQkJr6G3d4fqvqSkV2C1jsMtt3yDnYEKWmwnjwAYhokK8h7/iROncPHirut7HOgz8N5tr/FX\n7vM6g+7uXThxAjhxAmhtLQMAdgIRQrMDePzxx2G32zF69Gh88MEHqsesWrUK+/btw4033ojf/OY3\nyM3N1XpZhmF0iGT0Ozs/wyefWNDdLRn2CtlRSrOTEOS10lEAra02LFv2CrKzD3J4KAJodgDLly9H\naWkpli5dqrp/7969+Mtf/oKPP/4Yhw8fxlNPPYWWlhatl2UYRmf4h3nKAWyU7e0N8jrUPqV58o4I\nLl7chYYG7zs8ItCGZhXQnDlzkJKSEnT/nj17sGzZMgDArFmz4PF4cOHCBdVjZ8z4AWbMWIH8/ArY\nbOWw2xu1No9hmBhRXe24bvyBQONdCKBM5bV3OyHhSZV9SkehHBE0orXVgsceq2F7ESZRzwF0dnYi\nPT3dt52WloaOjg6MGTMm4Nj33hsDefzvnXdegdX6Jid/GEbHSGGfw4c7ZO8qjbf3tztq1GJkZ2fi\n0qULsFhWYvjwm5GYeA133z0dLS3rcfnyEN++np5efPLJk7Iwktxc9SWL3W7A4eDRQDjEJAmszFBb\nLJYgR3Lyh2GMRGDYR0Lqyff12K3WelRVrRjUb9hub8T27V7H4E0kS3uCjwbuusvBHcYBEnUHkJqa\nivb2dt92R0cHUlNTgxxdcf3/gwD+3W8PJ38YRn/4h33kRt/720xKWgSrdRxSU4ejtPT+Qf9mi4vn\n+v7G62zKrl9P3NGA0+mE0+mMzMkoArhcLsrOzlbdZ7fb6YEHHiAioubmZpo1a5bqcQAIoOv/Nshe\nEwENBDzrt52U9DBlZ6+mwsIyqqtriMTHYBhmANTVNVBhYRnl5W2g5OSlKr/Vcho5chnZbOUR/23W\n1TWQzVZOKSmLZNcsU7TB+89mK4/otfWKFjOu2QEsXryYxo0bR0OHDqW0tDSqqamhHTt20I4dO3zH\nrFy5kqxWK02fPp2OHj2q3hBAZuSVX2hZCGdAZLU+y06AYWJAXV0DWa3Pxt34+rdD2WH02omUlEWU\nl7fB9J1ELQ5AVzOBc3NLYLEMu578UWqIK66/VsrLvNhs61Ff/3xM2sowomKzlcPhkP/+lDN3Aav1\nWVRVDT7cM1i8+YEDOHLkY7jdtf20qQxVVTZThoS0zATWlQOQN0X6cvuSP9Iswgr4TyppBODAyJEd\nmDUrjXMDjCGwn7Zj9q2zkZyY7HsNAE1nm1A8uRieyx7faz0gV/r89a+/UextRErKLzB9+hQkJl5D\naWlBTH+D/oloQLROoqYyOpEYgkSCUE3xH+5xOIgxHnUf1ZG72+3bdne7qWR3CdV+UOt7XfJWCdWe\nqKU2dxutqFtB7m431X1UR23uNqr7qM7vb+XbUW970N+ffuLtUm4gL28DpaQo8xLef3l5G+Laxmih\nxYwbohx0cfFcVFXZYLOtR1bWZ0hKkiaNqE0V34Tt2w/EvI0Mo8R+2g7PZQ8AYPats1H25zKc8ZyB\n/bTdewABBz454D2GAFiA2//pdsx/bT7WzV6H5MRkZI/OxvzX5iN7dDbsp+044zmDsj+X+UYMnsue\nvvNFCXWlTx9W67MoLS2Iahv6o7h4Lurrn4fTWYG77kpX7G0EUI7333fxhDElEXREmhhMUyRvP3Lk\nMqGTP4y+cXe7fT15IqI2dxtlv5JNx7uO+953uV2ECpDL7fK9lva73C5aUbfCNyI43nWcsl/JpjZ3\nm+r5o0VenpoqL3pKH634j1jMHyXQYsYN6QAkCguVw1Hzf9mMvgkW6qk5WuMz4pLBlwy4y+2ikrdK\nqGR3ic/oy48jogDnUHOshkreKgm4ViRDQ5Lc019yqZ+wTyjU5aLGaPtgEdYB6EWSxohLqNi+tL3k\nj0sCevaSwXd3u/1yAO5ud8BIQRoBKJ3DkjeX+K4d6dFA/73ofzNExypw9GK+fICwDoCov+RPAwFl\nNHLkMg4HMVFBaXjd3W6/3rz0Whm+qT1R6zP4khNxd7up9oNan9GXEsDyv5M7B/l1in5X5DtG3rZw\nRwTqo+tySklZqsuwTzDUP0cZpaQsNY1NENoByPH/sjkcxMQGeShHGduXeumSMVc6C6WBVo4o5Cog\n6TrybeXoIlIjArP0nEXIB7ADuI4R5GqM8VEaaSJSje1HOk7fX35BGSoarPGXl3gYNeph0/x+zJ4P\nYAcgI7RCyHg9GEZ/hFL3yGP7asdGqw3StuSIjncdDzg+lOMJzKc1UELCE4aM+wfDLKMaJewAVBAh\n9sfED7nBDRbblx8b6YlbaqOQNncbFf9XcVhy0cDfi/c3M2qUV1JtpLh/MIJ/xocNLRtnB6CCCLE/\nJjaoGVt3t5tqjtaE1duOBsFGJdsPbx+QQzJr71jOwEY5xrML7ACCYPbYHxMb+lP6xGIyVn+EykuE\nkosaWesfDnLVoFnyHOwA+kGE3g0TXdQmbUU7zq+FUJPM5MbfDFr/cDGLXdBiOw1RC0grw4Yp1ycF\ngEacOHGKF6BnVJHX8QGA5MRkPDHzCUyomoB7b70X2wq3ITkx2bdv07xNaDrbFK/m+uG57EHZn8uw\nad4mZCRnYJttG3p6ezChagLWzl7ra7d/jZ+5AGwA1iMlZRlstvUxKekcT9TtApCYeC3GLYkfQjiA\nVasKYbXKC1g1IiHhNVy8uAsNDRVwODZi9er97AQYH1LxNskJnPGcwff/+H0cf/I4jp47GnB8cmKy\nbko3N51twqZ5m3yGHgCGJQxDzYIarLSvxBnPGQBAT8/1ZRUTPcAkO7xO4HlMnz4B9fXPm9r4A+p2\nISlpETo7/yZOpzCCIxFNRLspZoz9MdElmNJHjyGfYARLDre527yqmEQ3oWiF938BfweSXcjK+ldK\nSjJmQliL7RTGAcgxS+yPiSyhEql6UPqEg9pnqvlDLf3TKitNuusRGrJwqp/xN3vcPxjqElFjOEMt\ntlOIEJASjv0xaoQK+7z67qsBOQG9hHxCUTy52C8UZLc3ovKZ9/H5/30bHxe/jmuHn0aS5QlkZ/8Y\nNtt6bPrpbGDS3+LY4vjgC4cpuHx5SIxbEluEdAAc+2PUkJK5ZX8uw/sX3sf81+aj7tE6TB8z3fe+\n3AkYkepqB1o71wL3bgV+cRx46Gfo/toWpKYOR+1bT6MxYa9vsRmRELZTGMGRiCZi3RQzxP6Y6CCv\nvS/HKGGfUMye94x/zH9kG+GpbLqz6EnD5DWiQeAkMeOEw7TYTt0uCh8rbLZyOBziLCDN9CFfmF3i\njOcMVu5diZeLXsbWpq0BahqjM2PxIry3+1XgsuwzjXkfeCoHrtUuZCRnxK1t8cZub8T27Qdw+fIQ\nXLrUAeBrGDFiNIYN68WqVYW6VUWZflH4aMIJYXEJpZBR228GAnq6iW4a8cgMWvWf62JSw8gIqI8G\n9BsV0GI7hcwByBE29sf4xfzbPG1YYV+BukfrMD55vN9+vUzwChe7vRE2Wzny8ytQXe3AkiWpsNnW\nY/a8/4P0f/k2frno3/Hc4n8DLMCa/WvguezxTSYTMR/gP0HOS2vrJmzffiBOLYoewoeA7PZGrF69\n3+8LHzv2cYwbl2iI4R+jnTZPGyZUTTBlCETt+bZay1BVZQMm/c0vBOa57MEaxxrcm3Yvjp4/arrw\n10DJz69AQ0NFwPt5eRVwOgPfjzdabKe69kkgJMO+fft6X+zv/PlkvPfef/iOaW0t8zuWMS7KuL/n\nsgcbGzeiZkGNKWP+wXuzgTmu5MRklM8t9zlDM92HwSBSVED4EBDgNez19c/D6azAzTePRVfXf/jt\nN+vwT0TkWn+pxwsCHpr6kGmknnIGo2/3XPZga9NWuFa7/EpGyPfbT9uj0k49ESgT90YFPv30oulq\nhwk/AlAi6oQQUZDH/e+85U6AgG22vsJuUszfCJO8BsJAe7PyAnLJicn4efHPffMgxieP99tvdoSK\nCkQoEa0ZvTTFyFPCmYEjaf1dble8mxIVpDr/6vNcAvXtoVYY08uaB/FC7zZBi+3UHAKqr69HZmYm\nJk2ahBdffDFgv9PpxMiRI5Gbm4vc3Fxs3BioudcTPEvYfChLOyvj/mYK+QB9iV+HYyM+/PBVdHc/\niqSkRb5yD2plnpUlIwBgfPJ4vFz0ckAZadEwdVRAi+fp7e0lq9VKLpeLrly5Qjk5OXTy5Em/Y/77\nv/+bFixY0O+5NDYlovAsYXMh1/O7u91UsrvEp3k3o9Y/Uj1W+aIyZrtHg4FHAEE4cuQIJk6ciIyM\nDAwdOhSLFy/G7t271ZyMlsvEHCkpnJp6M7q7d/jt44Sw8ZDH/d889aZf3N8sWn85keixymP+H376\nIdbNXueXIBclIQyoJ4Wt1mdRWloQpxZFDk0OoLOzE+np6b7ttLQ0dHZ2+h1jsVhw6NAh5OTkoKio\nCCdPntRyyZhi6qGfYCQnJmPt7LUo2VOC8rxyv3CGUSp7DpRIyBjli8rMvnU2tjRtwbrZ69B0tkm4\nSWLFxXNRVWWDzbYeeXkVyM39AUaM8GDr1oOGDwtrUgFZLJZ+j5kxYwba29tx4403Yt++fXjwwQdx\n+vRp1WMrKip8r/Pz85Gfn6+leZoRSQ9sduQSRzPq/eWsWlWI1tYyxeSvZ1Faev+AzyF3iPIR1NrZ\na/3UQqJQXDwXxcVzVSfWxVoR5HQ64XQ6I3MyLbGn5uZmstlsvu3KykravHlzyL/JyMigixcvBryv\nsSlRIbAmSAMlJT1M2dmrqbCwjHMBOkWpaJHi/rUf1Pq2zRjTlpQ/eXkbKDe3hGbMWEF5eRvIZiuP\nyLNqduXUQNBjPkCL7dRkda9evUq33XYbuVwu6unpUU0Cd3V10VdffUVERIcPH6bx48erN0SHDoCI\nE8JGRGnga0/Umr7QWbQLmHFC2Isei0fGzQEQEe3du5cmT55MVquVKisriYhox44dtGPHDiIievnl\nlykrK4tycnLonnvuoebmZvWG6NQBSOjR8zPBEc1gRfP5lDvUuo/qqM3d5ndPzeZMQ6FHOxBXBxAp\n9O4A9Oj5mdCIFLKI5vMpD6lJzqDN3eZ7XwQHK6HHhWO02E6uBTRAOCFsLJRJX7NN9lISzedTPklM\nSghvadqCrNFZwiWETacIiqAj0oSOmqKKHj0/04fovdRYP58ija6CoZeFY7TYTuHXAxgM8iXjEhOv\n4e67x6G5+Rx6ehJ43YA4I5+41HS2Cdmjs7GlaYuvd+q57DFVkTc1lM9naWlBVJ5H6V6vnb0WK+0r\n8fPin/sW0ZH2m/1eA/pZTpaXhIwDevH+TB+iJX7lss9YyZJFXEYzGHrJC2qxnbqxukZzAHpUAzDi\nhCbi1QHhqqF96MUGaLGdnAQOEy4TEX9Eq/IpJ17r1nLV0D7MUCOIHUCYsCoo/oi2upccPXVARFNc\nSSgVQcFKbeuaCI5ENKGjpgwIVgXpAynmXHOsxvSzfeXoJfygjPmLlANQIx55GS22k1VAGpCrLi5d\n6gDwNYwYMZoVQTGmzdPmW8g8Izkj3s2JCWpFyazWZ2PeA7WftmP2rbORnJjsew3ApwISRREEBPtO\nylBVZYvqd8IqoDjDiqD4IbLyJxoF37Qg+mggXqMyLbZTN1bXyA5AL8Nx0RDN4BihoyGaQ5YTL1mo\nFtvJSeAIoKeEnJlRqn6azjb5FikBYMrVveTES/kzGKSFd0RTBAHGFIawA4gARvzijYhc9SNtb2na\n4rcyldlW95JjhI6GqIogwJiyUHYAEcCIX7wRka9M1eZpE64Qmd47GqKvI2zEQnGsAooQrAiKHSKq\nfgD9KH+CIVcESc5g3ex1OPHpCd/oTRSHHUtFEKuAdIQREnVGRsQko56VP6EQ8buSiKUwRIvt1I3V\nNYsDYEVQZOEyz8buUIhSm0lJLBVBWmwn5wAijBESdUZCnviVVD9S4tfsqh/AGMqfYIicENZ7vkaC\nHUCEMcoXbxTkid+s0Vl+Nf6l/WZV/QDG7VDIE8IZyRmYmzEXa/avCSjeZ9aksFGEIepPFxM2q1YV\norW1LCBRV1p6fxxbZWzk2nLXapcQSUQJo3Yoms42+Tlqm9WGA60HsP8v+7Eoe5GfgzAjUqJ3+/b1\nsgV69JGsl8MqoCgQq5WZREG+AtXWpq3CKEkA/St/BoPo32N1tSMqqwdqsZ3sAGJANL98M8JyQn/M\n1KEQUcIbbUkoy0B1jNFVHPFAru6p+6iO2txtATV/zFrmmSg+JYVjgaiy0GgrA7XYTt1YXbM6AJaF\nhoeoxsKsHQbRCvfJibYkVIvtZBVQlDGqiiPeiFpUzMiyz1DIk8KS8kcu4TWzIkjPiXx2AFFGz1++\nnhFVQ27WDoN8LWEpjyO9L+V55EX9zISuJaERGYNEAB01JaLw0pGDR+RwgSghQ9FCfHV1DWSzlUel\nhIcW28kqoBjAheJCI1f9SNvZo7Nx4tMTvkleoiwtaCbZZ3+IqAiSiKQyMK4qoH379tGUKVNo4sSJ\ntHnzZtVjSktLaeLEiTR9+nQ6duyY6jERaIruMWuCTysi9/jViGZvUS+INgKQE2k7oMV2arK6vb29\nZLVayeVy0ZUrVygnJ4dOnjzpd4zdbqcHHniAiIhaWlpo1qxZ6g25/iHMLPETZXgfDiIbBCLzSj/V\nUDr42hO1VPJWid93znZg4GhxAJqSwEeOHMHEiRORkZGBoUOHYvHixdi9e7ffMXv27MGyZcsAALNm\nzYLH48GFCxdUz2f2ZJBZE3yRQFTVD9AX9nE4NqKhoQIOx0asXr1flwuIRAK1MhGwAPv/sh8A24FY\noskBdHZ2Ij093bedlpaGzs7Ofo/p6OhQPZ/ZZ3iyIig4oqp+APNKP4MhVwQBXue/rXAbGs80CrHS\nm57sgKZicBaLZUDHkSJBEezvhr4zFC/9z0sAgPz8fOTn52tpnu7gQnF9hCr3IFX/NLMRkKOnHmG8\nEKngn1Y74HQ64XQ6I9IWTQ4gNTUV7e3tvu329nakpaWFPKajowOpqamq57s65yp+PO/Hpv3yjVIh\nMBbIa/rI6/xLRl963+yqH0BfPcJ4oRwBmtn5a7UDys7xc889F35jws4eENHVq1fptttuI5fLRT09\nPf0mgZubm0Mmgbnmi3kTf2qInviVEH2uiOgqMK12QIsZ16y93Lt3L02ePJmsVitVVlYSEdGOHTto\nx44dvmNWrlxJVquVpk+fTkePHlVvCCDUcn8sCfUi6pKBRMZd6zfSyJf9lF7LO39m7ghGwg7E1QFE\nCgBC9QhZEir2CIA7AOqINhqIhB3Q4gB0VQtIJCmgiIk/+2m7T90jT/x++OmHvsSvKOof0ZQ/A0W+\nBKgIiqB42wFdOQCRpIAiJv5EX+BdTrx/+HpGpI5g3O1A2GOHCIPrOQBRhn6iJv5EDvvI4RBgcER6\nRiJhB7SYcV0Vg3N3u331wqVZgJIU0IzFwJRL/d199zg0N58z/dKRIhcBkxCp6NtgkC8W33S2Cdmj\ns/3kwSLYgcEu+Wm6JSFFSwQRiZMUFKl31x8iFH0bLHJFkPSsiKIMJApPEqrFjOvSARCJZyhECAmI\n6NjliD7vIxxEsgPhdgJN6QCIxNKIR3vd0Hgg781J21JvTsLMGm85oozwooEodiDcTqAWB6ArFZAc\nkRRBgA7UAFFArvqRtiXVj0RyYrKp4rnBYNlneIhkB+KhDNOlA5AngjKSMzA3Yy7W7F/j9+WbbRFp\nXa8bGiaiabpDwbLPwaO0A2afKxKPTqAuHYCI9cKLi+eiqsoGm2098vIqYLOtN4UiRCRNdyjMOMKL\nNko7IM0dkc8VMVNHMC6dwLCDRxGmv6aIlAwiMk/CULTvLRiizvuIJCKICMJRhmkx44ZxAETiJIOM\nnDAUXcYXCpZ9ake0DsVAOoJCOACRvngjS0LlRl5S/YhU4luJWUZyeoI7gv7PkOkdgGhGxeiSUJGc\ndSiMPJLTKyI9WwPtCGpxALpMAiuRJ4MkKaGUDDJjQtjoCUNO/Hph6WdkkauCPvz0Q6ybvc5PFWSm\nhDAQG+WYIRyAfBFpSVq4pWkLskZnmVJaaHRJqEja7VCw9DOycEfQS0Q7glqHKZFisE0xexxQnjDM\nzS2h3NyndBtH5sSvOkbO5RgBs4eDBqoc02LGDekAzP7FyzFCHFm0HM1AYeln9BGpIxhMOSaUAxDN\n2BilFymSU+4PXus3Noj2zAVTlAnlAEQLNxhJEWT23thAMMKIzQyI1hEM9VwJ5QCUmL0XwCMAY2GU\n78voiNYRDPVcaXEAhlABhcLskkM9KoLki7sDXtXPGscazB0/V4iiXaFg5U9sEE0ZGK3nyvAOQC45\nXGlfiTOeMwH7jawNVisSt2RJGqqrHcjPr4DNVg67vTGmbVKWed7fuh8gwDbRBgDCLfAux+hzOIyK\n2TuCUXuuIjhK0UQ4TVEO9drcbZT9Sja1udtU95sBvcSYOeSjDit/4oPZn8dQz5UWM66rReEH2xRp\n8Xi5tz/jOYOVe1fi5aKXsbVpq+mGgjZbORyOjSrvr0d9/fMxbQsv7u7Fbm9EdbUDPT0JGDasF/fc\ncwtaWs6Hvcg3MzjkM4SlhePNGAZSLh5/993j0Nx8Dg7HpkHbTh8RclCaiWRTzKxGiZcqSLm8o7vb\nTSW7S6jmaI0pe1wDRS8jMpERcelR/+dO4CSwErOXIYhXjFke95eSviDgoakPCZ305Xo/8UeeEAYC\nlx41Y5kItecuHEzlAEQoFhUvVZB8ecc3T70JELDNtg3JiclCJ31Z9aM/RFiKNNhzN1jCPssXX3yB\nRYsW4cyZM8jIyMAbb7yB5OTAG5yRkYERI0ZgyJAhGDp0KI4cOaKpwaFQFosq+3OZr1iUtL1pnnav\nGU+kWPL27etx+fIQXLrUAeBr2Lr1IKqrHVi1qjBq8Wa50sK12uX3gxJlcXclrPrRJ6GeVTMQ7Lkb\nNOHGjtauXUsvvvgiERFt3ryZnnnmGdXjMjIy6OLFi/2eT0NTgiKmMiB68Wez389wYNWPPpGe1Zpj\nNVTyVklA/sro+YBI5QDC/sspU6ZQV1cXERGdP3+epkyZonpcRkYGff755/03JEr5aDMnhKM961S0\n2Zbhwks96gv5symJFSQnYKbnVnru4uIAkpOTfa+/+uorv205EyZMoDvuuIPuvPNO+uUvfxm8IVEc\nAZi1FxBtRZBo9VYGAy/1qF9EU6xpsZ0hcwAFBQXo6uoKeH/TJv84usVigcViUT1HU1MTxo0bh88+\n+wwFBQXIzMzEnDlzVI+tqKjwvc7Pz0d+fn6o5oVEnhAGgEPth7Bm/xpss20DAFPkA6Idf5Yn09bO\nXostTVv8kmmixv3t9kasXr3fT4XR2upNzLPeP/4on8nkxGSUzy03TT7A6XTC6XRG5mTheo4pU6bQ\n+fPniYjo3LlzQUNAcioqKmjbtm2q+zQ0RRURegGxij+bOYwWDlzwzViYPXelxXaGLQNduHAhdu7c\nCQDYuXMnHnzwwYBj/vGPf+Bvf/sbAODLL7+Ew+HAtGnTwr3koFBqg6VeQMn/KzFNrRBlnaDc3B9g\nxAgPtm49GFaNIGWRN8A7s/pHe39k2nkV4cDST+MggjRcE+F6josXL9K8efNo0qRJVFBQQG6316t2\ndnZSUVERERG1trZSTk4O5eTkUFZWFlVWVgY9n4amDAh5L6Dod0W+ekHy/UaOZ0dCESRibaVw4BGA\ncRBByKDFdpqyFIQSEQxbpIyS2R1luChX+Ro79ics/TQgZgwHsQPoB2U+gMjrBIr/q9g0D0IkFUEc\n8/dHbXQ1duzjvNSjQTHb863FdpqqFEQwlPkAABifPB4vF71smvrh4SqC1BZ32di4ETULajjmfx21\nuitdXTW4+eab4HRWoL7+eVb/GASzrx8yWIRwAGpID0LNwhpsbNgYYASN9hCEWyOIi7z1Dyd9zYE8\nIZyRnIHbgwJNAAAQgElEQVSfF/8c81+b73MCZiwa1x+GXg8gXJRzBCSjp5wjYLRRQbB64VKd+mB1\ngqT7cectd+LQ2UO+Im/SvqazTULq/SX0tAYDEz5mXT9Ek+2MUBhKM7FsirhzBIKrgswWF40kXO/H\n3Bj92ddiO4UcAahhttWtgvVaZyxajD//ZofpekHRQL7S16VLHbBYhmH48Jt5lS8TYYbRL48ANGJG\n6WMwVdDsec+YXhIbCXilL/NjlqJxWmyn8A7ArHMEQs0LMKPDizQ82cv8mCUUzA5AA2adI+DXg51U\nR0iso6Skhyk7ezUVFpZR9a4aQ8c9o0281l5m4osR8wFabKewMlCJ/uYIPDHzCb/9RpGIyusEZX79\n9xhSuA7d9CpOnHgJjsM/xNOHnkH17az1Dwav9CUeZpOGDwThHYAa0oNw/Mnj+P4fv29YnXBx8VzU\n1z+PW0en4ZqjCfhWGTDmfeDR+bj6m3dh/7WLtf4K7PZG2Gzl6Oz8DElJT/rti8Xay0x8kEvDH7r9\nIcACrNm/xjc/xki/+0ERwZGIJvTSlGA5geNdxw0VDpKHtnzhjDHHCRXw/i8LZ3DM30tg4rfBFzbj\ncg/mxsj5AC22Ux9Wl/TjANRyAse7vIbzeNdxv/f1bDjljqywsIwwso3wVLbX+BetICS6OaGpgBO/\njBwpH1BzrCbAJujpt6/FdnIISIEyJ+C57MGr776q+3CQsqZPcmIy1s1ehyVvLsGCH9yGof8yE3it\nDrgwHThYhCGFM3H208/DWjfArHDJB0ZCXjPo0NlDWONY47eGgJ5++5qIoCPShI6a4iNUOEhv8kll\nW6VtafSy/Y0astnKKSvrXykp6QlCoturDmJ9uw8eATBE6r+lkrdKqGR3iS6VgVpsp26srh4dQH/h\nIDWDG88HQ1nrXFrEXf7QspHzh+v8M0rUfvfubjfVHK3RZSiYHUCMCGZga47V+GYQyo+N9kMRjoOa\nPe8Z1rdfh+v8MwNFPqLW20RRdgAxoL8Qy5I3l8RkNBBqibuBhKhyFz3MI4Dr8GiIGQh6VwayA4gB\n/c0YlmKE0R4NBHsY3znzzoB6Juq93uWUm/sU5eVtoMLCMmF6vTzblxkIelcGsgOIA2ojgiV/XBIw\nGqg9UavJIQykVMVgH8a6ugay2cpDxL3FSArzCIAJB72Fg9gBxIFQE0fkigHptWSMpbxB7YlaX9XB\nuo/qfPuk19K+2g9q/cI8ytCTNAwNV50gohGUEr8+RRQnfZkBokdlIDuAOBNqNOByu1R7DMrys/Iq\npMH2ScZeciKR6IGIFgbh2b6MFvSiDJS3gx1AnBnINHJpVqG8x67UFstlm8p9ag+ZfEQhf+gG0+MI\nHAE0EFBGKSlLTZkPEHHEw0SPeCkD5c6FHYCOGMgkEsmYu9yugPKz8m01pxHpYaZ/j7iBAHMvgiLa\niIeJHgNVBkaioxYqF8gOQEcEGw3UflBLRP4xw4GMAGKRaJKSwikpi0zbO5bi/mb+jExsGagyUPkb\nlgtD5Pm+DQc3BM0LKiXffs6GHYB+Cabbr/2gdkA5gNoTtZp7DwPFrL3j/kc5nPhltBMqF6gW+pWM\nvfRbb3O3BbUJRIHzDyS7wA7AIMidgdzzh1IBSYY+FmoC9fh4A40a9bCh5wio5znKKSVlKSd+mYjR\nXy5QHvqV5w7UVIPByrmo5QLZATARQU0hk5CglEkaLydg1pENo18GMntYLd8XKi8YLBcYFwfwxhtv\n0NSpU+mGG26go0ePBj1u3759NGXKFJo4cSJt3rw5eEPYAegC+SSxUaOMWzZCXuTNyJ+DMSahSrYo\nhSEDGQGEygXGxQGcOnWKPvroI8rPzw/qAHp7e8lqtZLL5aIrV65QTk4OnTx5Ur0h7AB0R2DP2SsR\nHTlyma7DQQMbyXDcn4kNoYQhA80BhFISabGd6itgDIDMzMx+jzly5AgmTpyIjIwMAMDixYuxe/du\n3H777eFelokh/gujNwLYD2AT/vpXwOEAWlvLAHjXHtYT1dUOtLZukr0zF729wKhRi5GdnYnExGso\nLb1fd+1mzEnx5GK/7aazTdhWuA3Jicmwn7ZjW+E2AMBLzS/5XjedbQIA1D1ahxOfnvCdY9O8TWg6\n2+RbuEp57sEStgMYCJ2dnUhPT/dtp6Wl4fDhw9G8JBNBVq0qRGtr2XVj6gCwyW9/a+smbN++XjeG\n1G5vRHW1A4cPd6jsnYvs7INwOiti3SyG8UNutOWvK+6rUH1/fPJ43+tIGH05IR1AQUEBurq6At6v\nrKzEggUL+j25xWIJv2VM3JEM+/bt69HS0oG//lW+txGAAy0tHbDZyrFqVWFcHYHd3ojVq/dfd1bl\nqsckJl6LbaMYRueEdAAHDhzQdPLU1FS0t7f7ttvb25GWlhb0+IqKCt/r/Px85Ofna7o+o53i4rko\nLp4Lm60cDof0rv7CQf5hn0IAZZCPWKzWZ1Faen88msYwEcXpdMLpdEbmZFoTHPn5+fTuu++q7rt6\n9Srddttt5HK5qKenh5PABsY/saqfejqS2mfkyGWqWv+RI5ex1p8xNVpsZ9h/+eabb1JaWholJibS\nmDFj6P777ycios7OTioqKvIdt3fvXpo8eTJZrVaqrKwM3hB2ALpHkogGGluvwU1JWRTTCWN6dUoM\nE0u02E7L9RPEHYvFAp00hekHbzhoo+ydvpCQhNVahqoqW8RDQlKit6cnASdOnMLFi7tCtOFZVFWx\n2ocxN1psJzsAZtD4J1wBb9JV6RAcSElpx113pUcsQRx43Yrr/+TXPYCRI9tx993pKC0tYOPPmB52\nAEzMsdsbsX37AVy+PATvv++C273z+h5lT7wRSUmvwGodh1tu+UZYzkDq9f/P//wFbnetbI/S8Xix\n2dajvv75wX8ohjEg7ACYuOIfEpIb5cCwzNixJRg3bhhGjBiNYcN6fQ5BHtoZNqwX99xzC5qbz6Gz\n8zN88okF3d07oN7j57APIzZabGdUJ4IxYuA/YUz+SCknjzWiq2ssurr6RgfvvPMKRo9+FZ9+Ovy6\nkfe+f/Dga+jt3QF/hyKfmQwAXiPPM3wZJjzYATCakU8YO3LkY7jd0h7l4yV3CN7ee3f3Lpw5owzl\nOK4bf+U51PT99aiqWsFGn2HCgB0AExGkCWPeRK00GlD22IONDpSPoXxbfg7JyK9HSspZfPObt3KP\nn2E0wA6AiSjy0UBHx2f45JMnZaEduTEPZuSV28pe/9zrvf4SNvwMoxFOAjNRRa4WunSpA+fPj0BX\n138gdLK4EQkJr8nCQH1KotTU4SzvZBgZrAJiDIPkEDo6PpWpewClkb/77nFoaTmPy5eHXE/ustFn\nGDXYATCGRD46YCPPMOHBDoBhGEZQtNjOGyLcFoZhGMYgsANgGIYRFHYADMMwgsIOgGEYRlDYATAM\nwwgKOwCGYRhBYQfAMAwjKOwAGIZhBIUdAMMwjKCwA2AYhhEUdgAMwzCCwg6AYRhGUNgBMAzDCAo7\nAIZhGEFhB8AwDCMo7AAYhmEEhR0AwzCMoITtAH7/+98jKysLQ4YMwbFjx4Iel5GRgenTpyM3Nxff\n/OY3w70cwzAME2HCdgDTpk3Dn/70J8ydG3oNV4vFAqfTiffeew9HjhwJ93JC4XQ6490E3cD3og++\nF33wvYgMYTuAzMxMTJ48eUDH8lq/g4Mf7j74XvTB96IPvheRIeo5AIvFgm9/+9uYOXMmfvWrX0X7\ncgzDMMwASQi1s6CgAF1dXQHvV1ZWYsGCBQO6QFNTE8aNG4fPPvsMBQUFyMzMxJw5c8JrLcMwDBMx\nLKQxPnPffffhpz/9KWbMmNHvsc899xy+8Y1v4Omnnw5siMWipRkMwzDCEq4ZDzkC0Hrxf/zjH7h2\n7RqGDx+OL7/8Eg6HAxs2bBjUORiGYZjoEHYO4E9/+hPS09PR0tKC4uJiPPDAAwCAc+fOobi4GADQ\n1dWFOXPm4I477sCsWbMwf/58FBYWRqblDMMwjCY0h4AYhmEYYxL3mcD19fXIzMzEpEmT8OKLL8a7\nOTGlvb0d9913H7KyspCdnY3q6moAwBdffIGCggJMnjwZhYWF8Hg8cW5p7Lh27Rpyc3N9IgNR74XH\n48H3vvc93H777Zg6dSoOHz4s7L144YUXkJWVhWnTpuHRRx9FT0+PMPfi8ccfx5gxYzBt2jTfe6E+\n+wsvvIBJkyYhMzMTDoej3/PH1QFcu3YNP/rRj1BfX4+TJ0/i9ddfx6lTp+LZpJgydOhQ/OxnP8OH\nH36IlpYWvPLKKzh16hQ2b96MgoICnD59GvPmzcPmzZvj3dSYUVVVhalTp/pEAaLei9WrV6OoqAin\nTp3C+++/j8zMTCHvRVtbG371q1/h2LFj+OCDD3Dt2jXU1tYKcy+WL1+O+vp6v/eCffaTJ09i165d\nOHnyJOrr67FixQp89dVXoS9AceTQoUNks9l82y+88AK98MILcWxRfPnnf/5nOnDgAE2ZMoW6urqI\niOj8+fM0ZcqUOLcsNrS3t9O8efPo4MGDNH/+fCIiIe+Fx+OhCRMmBLwv4r24ePEiTZ48mb744gu6\nevUqzZ8/nxwOh1D3wuVyUXZ2tm872GevrKykzZs3+46z2WzU3Nwc8txxHQF0dnYiPT3dt52WlobO\nzs44tih+tLW14b333sOsWbNw4cIFjBkzBgAwZswYXLhwIc6tiw0/+clPsHXrVtxwQ99jKeK9cLlc\nuPnmm7F8+XLMmDEDP/zhD/Hll18KeS9uuukmPP3007j11ltxyy23IDk5GQUFBULeC4lgn/3cuXNI\nS0vzHTcQexpXB8Dafy9///vf8d3vfhdVVVUYPny43z6LxSLEfaqrq8Po0aORm5sbVBIsyr3o7e3F\nsWPHsGLFChw7dgxf//rXA0IcotyL1tZWvPTSS2hra8O5c+fw97//Hb/73e/8jhHlXqjR32fv777E\n1QGkpqaivb3dt93e3u7nwUTg6tWr+O53v4vHHnsMDz74IACvV5dmYJ8/fx6jR4+OZxNjwqFDh7Bn\nzx5MmDABjzzyCA4ePIjHHntMyHuRlpaGtLQ03HXXXQCA733vezh27BjGjh0r3L149913ce+992LU\nqFFISEjAQw89hObmZiHvhUSw34TSnnZ0dCA1NTXkueLqAGbOnImPP/4YbW1tuHLlCnbt2oWFCxfG\ns0kxhYhQUlKCqVOn4sc//rHv/YULF2Lnzp0AgJ07d/ocg5mprKxEe3s7XC4Xamtr8a1vfQu//e1v\nhbwXY8eORXp6Ok6fPg0AePvtt5GVlYUFCxYIdy8yMzPR0tKC7u5uEBHefvttTJ06Vch7IRHsN7Fw\n4ULU1tbiypUrcLlc+Pjjj/svwR/phMVg2bt3L02ePJmsVitVVlbGuzkx5Z133iGLxUI5OTl0xx13\n0B133EH79u2jixcv0rx582jSpElUUFBAbrc73k2NKU6nkxYsWEBEJOy9+N///V+aOXMmTZ8+nb7z\nne+Qx+MR9l68+OKLNHXqVMrOzqalS5fSlStXhLkXixcvpnHjxtHQoUMpLS2Nfv3rX4f87Js2bSKr\n1UpTpkyh+vr6fs/PE8EYhmEEJe4TwRiGYZj4wA6AYRhGUNgBMAzDCAo7AIZhGEFhB8AwDCMo7AAY\nhmEEhR0AwzCMoLADYBiGEZT/D72T+rREex4RAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xb1584c2c>"
       ]
      }
     ],
     "prompt_number": 71
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "> now, to compare how the different rk methods i.e 2,3,4 match up against one another."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rk4ylist = []\n",
      "rk4dylist = []\n",
      "rk4ylist.append(y)\n",
      "rk4dylist.append(dy)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 73
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "y = 1.\n",
      "dy = 0."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 74
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in xrange(iterations):\n",
      "    #k_update()\n",
      "    k1 = dy # y update\n",
      "    K1 = -y # dy update\n",
      "    k2 = dy + K1*h/2. # y update\n",
      "    K2 = -(y+h/2.)# dy update    \n",
      "    k3 = dy + K2*h/2. # y update\n",
      "    K3 = -(y+h/2.)# dy update\n",
      "    k4 = dy + K3*h # y update\n",
      "    K4 = -(y+h)# dy update\n",
      "    #print k1, k2, K1, K2\n",
      "    y += h*(k1+2*k2+2*k3+k4)/6.\n",
      "    dy += h*(K1+K2+K3+K4)/6.\n",
      "    #system_update()\n",
      "    #y += h*(k1+k2)/2.\n",
      "    #dy += h*(K1+K2)/2.\n",
      "    rk4ylist.append(y)\n",
      "    rk4dylist.append(dy)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 75
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(ylist,'o',dylist,'x')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 76,
       "text": [
        "[<matplotlib.lines.Line2D at 0xb1518b2c>,\n",
        " <matplotlib.lines.Line2D at 0xb151884c>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX90FPW5/98roSS2QCJXfphEgsuPSAIxiEXlQGJpspoA\n19oeQYtwMe1VoYH2CPi9JhziFSICvTVBW2xPesu3vRpsa4VvFsJiuZt4QgJX8KIIR2zcQBIIKu6W\n1oZA8Pn+scxmdnZ2k+zsj5n5PK9zOOzsTGY+Ozv7PJ/P87w/z8dCRASGYRhGOG6IdwMYhmGY+MAO\ngGEYRlDYATAMwwgKOwCGYRhBYQfAMAwjKOwAGIZhBEWzA3j88ccxZswYTJs2TXW/0+nEyJEjkZub\ni9zcXGzcuFHrJRmGYZgIkKD1BMuXL0dpaSmWLl0a9Ji8vDzs2bNH66UYhmGYCKJ5BDBnzhykpKSE\nPIbnmjEMw+iPqOcALBYLDh06hJycHBQVFeHkyZPRviTDMAwzADSHgPpjxowZaG9vx4033oh9+/bh\nwQcfxOnTp6N9WYZhGKYfou4Ahg8f7nv9wAMPYMWKFfjiiy9w0003+R1nsVii3RSGYRhTEm6YPeoh\noAsXLvgad+TIERBRgPGXICL+R4QNGzbEvQ16+cf3gu8F34vQ/7SgeQTwyCOPoKGhAZ9//jnS09Px\n3HPP4erVqwCAJ554An/4wx/wi1/8AgkJCbjxxhtRW1ur9ZIMwzBMBNDsAF5//fWQ+1euXImVK1dq\nvQzDMAwTYXgmsA7Jz8+PdxN0A9+LPvhe9MH3IjJYSGsQKUJYLBbN8SyGYRjR0GI7eQTAMAwjKOwA\nGIZhBIUdAMMwjKCwA2AYhhEUdgAMwzCCwg6AYRhGUNgBMAzDCAo7AIZhGEFhB8AwDCMo7AAYhmEE\nhR0AwzCMoLADYBiGERR2AAzDMILCDoBhGEZQ2AEwDMMICjsAhmEYQWEHwDAMIyjsABiGYQSFHQDD\nMIygsANgGIYRFHYADMMwgsIOgGEYRlDYATAMwwhKQrwbwDCMGNjtjaiudqCnJwGXLnUA+BpGjBiN\nYcN6cc89t6C5+ZzqvlWrClFcPDfezTclFiKieDcCACwWC3TSFIZhIozd3ojVq/ejtXUTgEYA+wFs\nur63EQkJr6G3d4fqvqSkV2C1jsMtt3yDnYEKWmwnjwAYhokK8h7/iROncPHirut7HOgz8N5tr/FX\n7vM6g+7uXThxAjhxAmhtLQMAdgIRQrMDePzxx2G32zF69Gh88MEHqsesWrUK+/btw4033ojf/OY3\nyM3N1XpZhmF0iGT0Ozs/wyefWNDdLRn2CtlRSrOTEOS10lEAra02LFv2CrKzD3J4KAJodgDLly9H\naWkpli5dqrp/7969+Mtf/oKPP/4Yhw8fxlNPPYWWlhatl2UYRmf4h3nKAWyU7e0N8jrUPqV58o4I\nLl7chYYG7zs8ItCGZhXQnDlzkJKSEnT/nj17sGzZMgDArFmz4PF4cOHCBdVjZ8z4AWbMWIH8/ArY\nbOWw2xu1No9hmBhRXe24bvyBQONdCKBM5bV3OyHhSZV9SkehHBE0orXVgsceq2F7ESZRzwF0dnYi\nPT3dt52WloaOjg6MGTMm4Nj33hsDefzvnXdegdX6Jid/GEbHSGGfw4c7ZO8qjbf3tztq1GJkZ2fi\n0qULsFhWYvjwm5GYeA133z0dLS3rcfnyEN++np5efPLJk7Iwktxc9SWL3W7A4eDRQDjEJAmszFBb\nLJYgR3Lyh2GMRGDYR0Lqyff12K3WelRVrRjUb9hub8T27V7H4E0kS3uCjwbuusvBHcYBEnUHkJqa\nivb2dt92R0cHUlNTgxxdcf3/gwD+3W8PJ38YRn/4h33kRt/720xKWgSrdRxSU4ejtPT+Qf9mi4vn\n+v7G62zKrl9P3NGA0+mE0+mMzMkoArhcLsrOzlbdZ7fb6YEHHiAioubmZpo1a5bqcQAIoOv/Nshe\nEwENBDzrt52U9DBlZ6+mwsIyqqtriMTHYBhmANTVNVBhYRnl5W2g5OSlKr/Vcho5chnZbOUR/23W\n1TWQzVZOKSmLZNcsU7TB+89mK4/otfWKFjOu2QEsXryYxo0bR0OHDqW0tDSqqamhHTt20I4dO3zH\nrFy5kqxWK02fPp2OHj2q3hBAZuSVX2hZCGdAZLU+y06AYWJAXV0DWa3Pxt34+rdD2WH02omUlEWU\nl7fB9J1ELQ5AVzOBc3NLYLEMu578UWqIK66/VsrLvNhs61Ff/3xM2sowomKzlcPhkP/+lDN3Aav1\nWVRVDT7cM1i8+YEDOHLkY7jdtf20qQxVVTZThoS0zATWlQOQN0X6cvuSP9Iswgr4TyppBODAyJEd\nmDUrjXMDjCGwn7Zj9q2zkZyY7HsNAE1nm1A8uRieyx7faz0gV/r89a+/UextRErKLzB9+hQkJl5D\naWlBTH+D/oloQLROoqYyOpEYgkSCUE3xH+5xOIgxHnUf1ZG72+3bdne7qWR3CdV+UOt7XfJWCdWe\nqKU2dxutqFtB7m431X1UR23uNqr7qM7vb+XbUW970N+ffuLtUm4gL28DpaQo8xLef3l5G+Laxmih\nxYwbohx0cfFcVFXZYLOtR1bWZ0hKkiaNqE0V34Tt2w/EvI0Mo8R+2g7PZQ8AYPats1H25zKc8ZyB\n/bTdewABBz454D2GAFiA2//pdsx/bT7WzV6H5MRkZI/OxvzX5iN7dDbsp+044zmDsj+X+UYMnsue\nvvNFCXWlTx9W67MoLS2Iahv6o7h4Lurrn4fTWYG77kpX7G0EUI7333fxhDElEXREmhhMUyRvP3Lk\nMqGTP4y+cXe7fT15IqI2dxtlv5JNx7uO+953uV2ECpDL7fK9lva73C5aUbfCNyI43nWcsl/JpjZ3\nm+r5o0VenpoqL3pKH634j1jMHyXQYsYN6QAkCguVw1Hzf9mMvgkW6qk5WuMz4pLBlwy4y+2ikrdK\nqGR3ic/oy48jogDnUHOshkreKgm4ViRDQ5Lc019yqZ+wTyjU5aLGaPtgEdYB6EWSxohLqNi+tL3k\nj0sCevaSwXd3u/1yAO5ud8BIQRoBKJ3DkjeX+K4d6dFA/73ofzNExypw9GK+fICwDoCov+RPAwFl\nNHLkMg4HMVFBaXjd3W6/3rz0Whm+qT1R6zP4khNxd7up9oNan9GXEsDyv5M7B/l1in5X5DtG3rZw\nRwTqo+tySklZqsuwTzDUP0cZpaQsNY1NENoByPH/sjkcxMQGeShHGduXeumSMVc6C6WBVo4o5Cog\n6TrybeXoIlIjArP0nEXIB7ADuI4R5GqM8VEaaSJSje1HOk7fX35BGSoarPGXl3gYNeph0/x+zJ4P\nYAcgI7RCyHg9GEZ/hFL3yGP7asdGqw3StuSIjncdDzg+lOMJzKc1UELCE4aM+wfDLKMaJewAVBAh\n9sfED7nBDRbblx8b6YlbaqOQNncbFf9XcVhy0cDfi/c3M2qUV1JtpLh/MIJ/xocNLRtnB6CCCLE/\nJjaoGVt3t5tqjtaE1duOBsFGJdsPbx+QQzJr71jOwEY5xrML7ACCYPbYHxMb+lP6xGIyVn+EykuE\nkosaWesfDnLVoFnyHOwA+kGE3g0TXdQmbUU7zq+FUJPM5MbfDFr/cDGLXdBiOw1RC0grw4Yp1ycF\ngEacOHGKF6BnVJHX8QGA5MRkPDHzCUyomoB7b70X2wq3ITkx2bdv07xNaDrbFK/m+uG57EHZn8uw\nad4mZCRnYJttG3p6ezChagLWzl7ra7d/jZ+5AGwA1iMlZRlstvUxKekcT9TtApCYeC3GLYkfQjiA\nVasKYbXKC1g1IiHhNVy8uAsNDRVwODZi9er97AQYH1LxNskJnPGcwff/+H0cf/I4jp47GnB8cmKy\nbko3N51twqZ5m3yGHgCGJQxDzYIarLSvxBnPGQBAT8/1ZRUTPcAkO7xO4HlMnz4B9fXPm9r4A+p2\nISlpETo7/yZOpzCCIxFNRLspZoz9MdElmNJHjyGfYARLDre527yqmEQ3oWiF938BfweSXcjK+ldK\nSjJmQliL7RTGAcgxS+yPiSyhEql6UPqEg9pnqvlDLf3TKitNuusRGrJwqp/xN3vcPxjqElFjOEMt\ntlOIEJASjv0xaoQK+7z67qsBOQG9hHxCUTy52C8UZLc3ovKZ9/H5/30bHxe/jmuHn0aS5QlkZ/8Y\nNtt6bPrpbGDS3+LY4vjgC4cpuHx5SIxbEluEdAAc+2PUkJK5ZX8uw/sX3sf81+aj7tE6TB8z3fe+\n3AkYkepqB1o71wL3bgV+cRx46Gfo/toWpKYOR+1bT6MxYa9vsRmRELZTGMGRiCZi3RQzxP6Y6CCv\nvS/HKGGfUMye94x/zH9kG+GpbLqz6EnD5DWiQeAkMeOEw7TYTt0uCh8rbLZyOBziLCDN9CFfmF3i\njOcMVu5diZeLXsbWpq0BahqjM2PxIry3+1XgsuwzjXkfeCoHrtUuZCRnxK1t8cZub8T27Qdw+fIQ\nXLrUAeBrGDFiNIYN68WqVYW6VUWZflH4aMIJYXEJpZBR228GAnq6iW4a8cgMWvWf62JSw8gIqI8G\n9BsV0GI7hcwByBE29sf4xfzbPG1YYV+BukfrMD55vN9+vUzwChe7vRE2Wzny8ytQXe3AkiWpsNnW\nY/a8/4P0f/k2frno3/Hc4n8DLMCa/WvguezxTSYTMR/gP0HOS2vrJmzffiBOLYoewoeA7PZGrF69\n3+8LHzv2cYwbl2iI4R+jnTZPGyZUTTBlCETt+bZay1BVZQMm/c0vBOa57MEaxxrcm3Yvjp4/arrw\n10DJz69AQ0NFwPt5eRVwOgPfjzdabKe69kkgJMO+fft6X+zv/PlkvPfef/iOaW0t8zuWMS7KuL/n\nsgcbGzeiZkGNKWP+wXuzgTmu5MRklM8t9zlDM92HwSBSVED4EBDgNez19c/D6azAzTePRVfXf/jt\nN+vwT0TkWn+pxwsCHpr6kGmknnIGo2/3XPZga9NWuFa7/EpGyPfbT9uj0k49ESgT90YFPv30oulq\nhwk/AlAi6oQQUZDH/e+85U6AgG22vsJuUszfCJO8BsJAe7PyAnLJicn4efHPffMgxieP99tvdoSK\nCkQoEa0ZvTTFyFPCmYEjaf1dble8mxIVpDr/6vNcAvXtoVYY08uaB/FC7zZBi+3UHAKqr69HZmYm\nJk2ahBdffDFgv9PpxMiRI5Gbm4vc3Fxs3BioudcTPEvYfChLOyvj/mYK+QB9iV+HYyM+/PBVdHc/\niqSkRb5yD2plnpUlIwBgfPJ4vFz0ckAZadEwdVRAi+fp7e0lq9VKLpeLrly5Qjk5OXTy5Em/Y/77\nv/+bFixY0O+5NDYlovAsYXMh1/O7u91UsrvEp3k3o9Y/Uj1W+aIyZrtHg4FHAEE4cuQIJk6ciIyM\nDAwdOhSLFy/G7t271ZyMlsvEHCkpnJp6M7q7d/jt44Sw8ZDH/d889aZf3N8sWn85keixymP+H376\nIdbNXueXIBclIQyoJ4Wt1mdRWloQpxZFDk0OoLOzE+np6b7ttLQ0dHZ2+h1jsVhw6NAh5OTkoKio\nCCdPntRyyZhi6qGfYCQnJmPt7LUo2VOC8rxyv3CGUSp7DpRIyBjli8rMvnU2tjRtwbrZ69B0tkm4\nSWLFxXNRVWWDzbYeeXkVyM39AUaM8GDr1oOGDwtrUgFZLJZ+j5kxYwba29tx4403Yt++fXjwwQdx\n+vRp1WMrKip8r/Pz85Gfn6+leZoRSQ9sduQSRzPq/eWsWlWI1tYyxeSvZ1Faev+AzyF3iPIR1NrZ\na/3UQqJQXDwXxcVzVSfWxVoR5HQ64XQ6I3MyLbGn5uZmstlsvu3KykravHlzyL/JyMigixcvBryv\nsSlRIbAmSAMlJT1M2dmrqbCwjHMBOkWpaJHi/rUf1Pq2zRjTlpQ/eXkbKDe3hGbMWEF5eRvIZiuP\nyLNqduXUQNBjPkCL7dRkda9evUq33XYbuVwu6unpUU0Cd3V10VdffUVERIcPH6bx48erN0SHDoCI\nE8JGRGnga0/Umr7QWbQLmHFC2Isei0fGzQEQEe3du5cmT55MVquVKisriYhox44dtGPHDiIievnl\nlykrK4tycnLonnvuoebmZvWG6NQBSOjR8zPBEc1gRfP5lDvUuo/qqM3d5ndPzeZMQ6FHOxBXBxAp\n9O4A9Oj5mdCIFLKI5vMpD6lJzqDN3eZ7XwQHK6HHhWO02E6uBTRAOCFsLJRJX7NN9lISzedTPklM\nSghvadqCrNFZwiWETacIiqAj0oSOmqKKHj0/04fovdRYP58ija6CoZeFY7TYTuHXAxgM8iXjEhOv\n4e67x6G5+Rx6ehJ43YA4I5+41HS2Cdmjs7GlaYuvd+q57DFVkTc1lM9naWlBVJ5H6V6vnb0WK+0r\n8fPin/sW0ZH2m/1eA/pZTpaXhIwDevH+TB+iJX7lss9YyZJFXEYzGHrJC2qxnbqxukZzAHpUAzDi\nhCbi1QHhqqF96MUGaLGdnAQOEy4TEX9Eq/IpJ17r1nLV0D7MUCOIHUCYsCoo/oi2upccPXVARFNc\nSSgVQcFKbeuaCI5ENKGjpgwIVgXpAynmXHOsxvSzfeXoJfygjPmLlANQIx55GS22k1VAGpCrLi5d\n6gDwNYwYMZoVQTGmzdPmW8g8Izkj3s2JCWpFyazWZ2PeA7WftmP2rbORnJjsew3ApwISRREEBPtO\nylBVZYvqd8IqoDjDiqD4IbLyJxoF37Qg+mggXqMyLbZTN1bXyA5AL8Nx0RDN4BihoyGaQ5YTL1mo\nFtvJSeAIoKeEnJlRqn6azjb5FikBYMrVveTES/kzGKSFd0RTBAHGFIawA4gARvzijYhc9SNtb2na\n4rcyldlW95JjhI6GqIogwJiyUHYAEcCIX7wRka9M1eZpE64Qmd47GqKvI2zEQnGsAooQrAiKHSKq\nfgD9KH+CIVcESc5g3ex1OPHpCd/oTRSHHUtFEKuAdIQREnVGRsQko56VP6EQ8buSiKUwRIvt1I3V\nNYsDYEVQZOEyz8buUIhSm0lJLBVBWmwn5wAijBESdUZCnviVVD9S4tfsqh/AGMqfYIicENZ7vkaC\nHUCEMcoXbxTkid+s0Vl+Nf6l/WZV/QDG7VDIE8IZyRmYmzEXa/avCSjeZ9aksFGEIepPFxM2q1YV\norW1LCBRV1p6fxxbZWzk2nLXapcQSUQJo3Yoms42+Tlqm9WGA60HsP8v+7Eoe5GfgzAjUqJ3+/b1\nsgV69JGsl8MqoCgQq5WZREG+AtXWpq3CKEkA/St/BoPo32N1tSMqqwdqsZ3sAGJANL98M8JyQn/M\n1KEQUcIbbUkoy0B1jNFVHPFAru6p+6iO2txtATV/zFrmmSg+JYVjgaiy0GgrA7XYTt1YXbM6AJaF\nhoeoxsKsHQbRCvfJibYkVIvtZBVQlDGqiiPeiFpUzMiyz1DIk8KS8kcu4TWzIkjPiXx2AFFGz1++\nnhFVQ27WDoN8LWEpjyO9L+V55EX9zISuJaERGYNEAB01JaLw0pGDR+RwgSghQ9FCfHV1DWSzlUel\nhIcW28kqoBjAheJCI1f9SNvZo7Nx4tMTvkleoiwtaCbZZ3+IqAiSiKQyMK4qoH379tGUKVNo4sSJ\ntHnzZtVjSktLaeLEiTR9+nQ6duyY6jERaIruMWuCTysi9/jViGZvUS+INgKQE2k7oMV2arK6vb29\nZLVayeVy0ZUrVygnJ4dOnjzpd4zdbqcHHniAiIhaWlpo1qxZ6g25/iHMLPETZXgfDiIbBCLzSj/V\nUDr42hO1VPJWid93znZg4GhxAJqSwEeOHMHEiRORkZGBoUOHYvHixdi9e7ffMXv27MGyZcsAALNm\nzYLH48GFCxdUz2f2ZJBZE3yRQFTVD9AX9nE4NqKhoQIOx0asXr1flwuIRAK1MhGwAPv/sh8A24FY\noskBdHZ2Ij093bedlpaGzs7Ofo/p6OhQPZ/ZZ3iyIig4oqp+APNKP4MhVwQBXue/rXAbGs80CrHS\nm57sgKZicBaLZUDHkSJBEezvhr4zFC/9z0sAgPz8fOTn52tpnu7gQnF9hCr3IFX/NLMRkKOnHmG8\nEKngn1Y74HQ64XQ6I9IWTQ4gNTUV7e3tvu329nakpaWFPKajowOpqamq57s65yp+PO/Hpv3yjVIh\nMBbIa/rI6/xLRl963+yqH0BfPcJ4oRwBmtn5a7UDys7xc889F35jws4eENHVq1fptttuI5fLRT09\nPf0mgZubm0Mmgbnmi3kTf2qInviVEH2uiOgqMK12QIsZ16y93Lt3L02ePJmsVitVVlYSEdGOHTto\nx44dvmNWrlxJVquVpk+fTkePHlVvCCDUcn8sCfUi6pKBRMZd6zfSyJf9lF7LO39m7ghGwg7E1QFE\nCgBC9QhZEir2CIA7AOqINhqIhB3Q4gB0VQtIJCmgiIk/+2m7T90jT/x++OmHvsSvKOof0ZQ/A0W+\nBKgIiqB42wFdOQCRpIAiJv5EX+BdTrx/+HpGpI5g3O1A2GOHCIPrOQBRhn6iJv5EDvvI4RBgcER6\nRiJhB7SYcV0Vg3N3u331wqVZgJIU0IzFwJRL/d199zg0N58z/dKRIhcBkxCp6NtgkC8W33S2Cdmj\ns/3kwSLYgcEu+Wm6JSFFSwQRiZMUFKl31x8iFH0bLHJFkPSsiKIMJApPEqrFjOvSARCJZyhECAmI\n6NjliD7vIxxEsgPhdgJN6QCIxNKIR3vd0Hgg781J21JvTsLMGm85oozwooEodiDcTqAWB6ArFZAc\nkRRBgA7UAFFArvqRtiXVj0RyYrKp4rnBYNlneIhkB+KhDNOlA5AngjKSMzA3Yy7W7F/j9+WbbRFp\nXa8bGiaiabpDwbLPwaO0A2afKxKPTqAuHYCI9cKLi+eiqsoGm2098vIqYLOtN4UiRCRNdyjMOMKL\nNko7IM0dkc8VMVNHMC6dwLCDRxGmv6aIlAwiMk/CULTvLRiizvuIJCKICMJRhmkx44ZxAETiJIOM\nnDAUXcYXCpZ9ake0DsVAOoJCOACRvngjS0LlRl5S/YhU4luJWUZyeoI7gv7PkOkdgGhGxeiSUJGc\ndSiMPJLTKyI9WwPtCGpxALpMAiuRJ4MkKaGUDDJjQtjoCUNO/Hph6WdkkauCPvz0Q6ybvc5PFWSm\nhDAQG+WYIRyAfBFpSVq4pWkLskZnmVJaaHRJqEja7VCw9DOycEfQS0Q7glqHKZFisE0xexxQnjDM\nzS2h3NyndBtH5sSvOkbO5RgBs4eDBqoc02LGDekAzP7FyzFCHFm0HM1AYeln9BGpIxhMOSaUAxDN\n2BilFymSU+4PXus3Noj2zAVTlAnlAEQLNxhJEWT23thAMMKIzQyI1hEM9VwJ5QCUmL0XwCMAY2GU\n78voiNYRDPVcaXEAhlABhcLskkM9KoLki7sDXtXPGscazB0/V4iiXaFg5U9sEE0ZGK3nyvAOQC45\nXGlfiTOeMwH7jawNVisSt2RJGqqrHcjPr4DNVg67vTGmbVKWed7fuh8gwDbRBgDCLfAux+hzOIyK\n2TuCUXuuIjhK0UQ4TVEO9drcbZT9Sja1udtU95sBvcSYOeSjDit/4oPZn8dQz5UWM66rReEH2xRp\n8Xi5tz/jOYOVe1fi5aKXsbVpq+mGgjZbORyOjSrvr0d9/fMxbQsv7u7Fbm9EdbUDPT0JGDasF/fc\ncwtaWs6Hvcg3MzjkM4SlhePNGAZSLh5/993j0Nx8Dg7HpkHbTh8RclCaiWRTzKxGiZcqSLm8o7vb\nTSW7S6jmaI0pe1wDRS8jMpERcelR/+dO4CSwErOXIYhXjFke95eSviDgoakPCZ305Xo/8UeeEAYC\nlx41Y5kItecuHEzlAEQoFhUvVZB8ecc3T70JELDNtg3JiclCJ31Z9aM/RFiKNNhzN1jCPssXX3yB\nRYsW4cyZM8jIyMAbb7yB5OTAG5yRkYERI0ZgyJAhGDp0KI4cOaKpwaFQFosq+3OZr1iUtL1pnnav\nGU+kWPL27etx+fIQXLrUAeBr2Lr1IKqrHVi1qjBq8Wa50sK12uX3gxJlcXclrPrRJ6GeVTMQ7Lkb\nNOHGjtauXUsvvvgiERFt3ryZnnnmGdXjMjIy6OLFi/2eT0NTgiKmMiB68Wez389wYNWPPpGe1Zpj\nNVTyVklA/sro+YBI5QDC/sspU6ZQV1cXERGdP3+epkyZonpcRkYGff755/03JEr5aDMnhKM961S0\n2Zbhwks96gv5symJFSQnYKbnVnru4uIAkpOTfa+/+uorv205EyZMoDvuuIPuvPNO+uUvfxm8IVEc\nAZi1FxBtRZBo9VYGAy/1qF9EU6xpsZ0hcwAFBQXo6uoKeH/TJv84usVigcViUT1HU1MTxo0bh88+\n+wwFBQXIzMzEnDlzVI+tqKjwvc7Pz0d+fn6o5oVEnhAGgEPth7Bm/xpss20DAFPkA6Idf5Yn09bO\nXostTVv8kmmixv3t9kasXr3fT4XR2upNzLPeP/4on8nkxGSUzy03TT7A6XTC6XRG5mTheo4pU6bQ\n+fPniYjo3LlzQUNAcioqKmjbtm2q+zQ0RRURegGxij+bOYwWDlzwzViYPXelxXaGLQNduHAhdu7c\nCQDYuXMnHnzwwYBj/vGPf+Bvf/sbAODLL7+Ew+HAtGnTwr3koFBqg6VeQMn/KzFNrRBlnaDc3B9g\nxAgPtm49GFaNIGWRN8A7s/pHe39k2nkV4cDST+MggjRcE+F6josXL9K8efNo0qRJVFBQQG6316t2\ndnZSUVERERG1trZSTk4O5eTkUFZWFlVWVgY9n4amDAh5L6Dod0W+ekHy/UaOZ0dCESRibaVw4BGA\ncRBByKDFdpqyFIQSEQxbpIyS2R1luChX+Ro79ics/TQgZgwHsQPoB2U+gMjrBIr/q9g0D0IkFUEc\n8/dHbXQ1duzjvNSjQTHb863FdpqqFEQwlPkAABifPB4vF71smvrh4SqC1BZ32di4ETULajjmfx21\nuitdXTW4+eab4HRWoL7+eVb/GASzrx8yWIRwAGpID0LNwhpsbNgYYASN9hCEWyOIi7z1Dyd9zYE8\nIZyRnIHbgwJNAAAQgElEQVSfF/8c81+b73MCZiwa1x+GXg8gXJRzBCSjp5wjYLRRQbB64VKd+mB1\ngqT7cectd+LQ2UO+Im/SvqazTULq/SX0tAYDEz5mXT9Ek+2MUBhKM7FsirhzBIKrgswWF40kXO/H\n3Bj92ddiO4UcAahhttWtgvVaZyxajD//ZofpekHRQL7S16VLHbBYhmH48Jt5lS8TYYbRL48ANGJG\n6WMwVdDsec+YXhIbCXilL/NjlqJxWmyn8A7ArHMEQs0LMKPDizQ82cv8mCUUzA5AA2adI+DXg51U\nR0iso6Skhyk7ezUVFpZR9a4aQ8c9o0281l5m4osR8wFabKewMlCJ/uYIPDHzCb/9RpGIyusEZX79\n9xhSuA7d9CpOnHgJjsM/xNOHnkH17az1Dwav9CUeZpOGDwThHYAa0oNw/Mnj+P4fv29YnXBx8VzU\n1z+PW0en4ZqjCfhWGTDmfeDR+bj6m3dh/7WLtf4K7PZG2Gzl6Oz8DElJT/rti8Xay0x8kEvDH7r9\nIcACrNm/xjc/xki/+0ERwZGIJvTSlGA5geNdxw0VDpKHtnzhjDHHCRXw/i8LZ3DM30tg4rfBFzbj\ncg/mxsj5AC22Ux9Wl/TjANRyAse7vIbzeNdxv/f1bDjljqywsIwwso3wVLbX+BetICS6OaGpgBO/\njBwpH1BzrCbAJujpt6/FdnIISIEyJ+C57MGr776q+3CQsqZPcmIy1s1ehyVvLsGCH9yGof8yE3it\nDrgwHThYhCGFM3H208/DWjfArHDJB0ZCXjPo0NlDWONY47eGgJ5++5qIoCPShI6a4iNUOEhv8kll\nW6VtafSy/Y0astnKKSvrXykp6QlCoturDmJ9uw8eATBE6r+lkrdKqGR3iS6VgVpsp26srh4dQH/h\nIDWDG88HQ1nrXFrEXf7QspHzh+v8M0rUfvfubjfVHK3RZSiYHUCMCGZga47V+GYQyo+N9kMRjoOa\nPe8Z1rdfh+v8MwNFPqLW20RRdgAxoL8Qy5I3l8RkNBBqibuBhKhyFz3MI4Dr8GiIGQh6VwayA4gB\n/c0YlmKE0R4NBHsY3znzzoB6Juq93uWUm/sU5eVtoMLCMmF6vTzblxkIelcGsgOIA2ojgiV/XBIw\nGqg9UavJIQykVMVgH8a6ugay2cpDxL3FSArzCIAJB72Fg9gBxIFQE0fkigHptWSMpbxB7YlaX9XB\nuo/qfPuk19K+2g9q/cI8ytCTNAwNV50gohGUEr8+RRQnfZkBokdlIDuAOBNqNOByu1R7DMrys/Iq\npMH2ScZeciKR6IGIFgbh2b6MFvSiDJS3gx1AnBnINHJpVqG8x67UFstlm8p9ag+ZfEQhf+gG0+MI\nHAE0EFBGKSlLTZkPEHHEw0SPeCkD5c6FHYCOGMgkEsmYu9yugPKz8m01pxHpYaZ/j7iBAHMvgiLa\niIeJHgNVBkaioxYqF8gOQEcEGw3UflBLRP4xw4GMAGKRaJKSwikpi0zbO5bi/mb+jExsGagyUPkb\nlgtD5Pm+DQc3BM0LKiXffs6GHYB+Cabbr/2gdkA5gNoTtZp7DwPFrL3j/kc5nPhltBMqF6gW+pWM\nvfRbb3O3BbUJRIHzDyS7wA7AIMidgdzzh1IBSYY+FmoC9fh4A40a9bCh5wio5znKKSVlKSd+mYjR\nXy5QHvqV5w7UVIPByrmo5QLZATARQU0hk5CglEkaLydg1pENo18GMntYLd8XKi8YLBcYFwfwxhtv\n0NSpU+mGG26go0ePBj1u3759NGXKFJo4cSJt3rw5eEPYAegC+SSxUaOMWzZCXuTNyJ+DMSahSrYo\nhSEDGQGEygXGxQGcOnWKPvroI8rPzw/qAHp7e8lqtZLL5aIrV65QTk4OnTx5Ur0h7AB0R2DP2SsR\nHTlyma7DQQMbyXDcn4kNoYQhA80BhFISabGd6itgDIDMzMx+jzly5AgmTpyIjIwMAMDixYuxe/du\n3H777eFelokh/gujNwLYD2AT/vpXwOEAWlvLAHjXHtYT1dUOtLZukr0zF729wKhRi5GdnYnExGso\nLb1fd+1mzEnx5GK/7aazTdhWuA3Jicmwn7ZjW+E2AMBLzS/5XjedbQIA1D1ahxOfnvCdY9O8TWg6\n2+RbuEp57sEStgMYCJ2dnUhPT/dtp6Wl4fDhw9G8JBNBVq0qRGtr2XVj6gCwyW9/a+smbN++XjeG\n1G5vRHW1A4cPd6jsnYvs7INwOiti3SyG8UNutOWvK+6rUH1/fPJ43+tIGH05IR1AQUEBurq6At6v\nrKzEggUL+j25xWIJv2VM3JEM+/bt69HS0oG//lW+txGAAy0tHbDZyrFqVWFcHYHd3ojVq/dfd1bl\nqsckJl6LbaMYRueEdAAHDhzQdPLU1FS0t7f7ttvb25GWlhb0+IqKCt/r/Px85Ofna7o+o53i4rko\nLp4Lm60cDof0rv7CQf5hn0IAZZCPWKzWZ1Faen88msYwEcXpdMLpdEbmZFoTHPn5+fTuu++q7rt6\n9Srddttt5HK5qKenh5PABsY/saqfejqS2mfkyGWqWv+RI5ex1p8xNVpsZ9h/+eabb1JaWholJibS\nmDFj6P777ycios7OTioqKvIdt3fvXpo8eTJZrVaqrKwM3hB2ALpHkogGGluvwU1JWRTTCWN6dUoM\nE0u02E7L9RPEHYvFAp00hekHbzhoo+ydvpCQhNVahqoqW8RDQlKit6cnASdOnMLFi7tCtOFZVFWx\n2ocxN1psJzsAZtD4J1wBb9JV6RAcSElpx113pUcsQRx43Yrr/+TXPYCRI9tx993pKC0tYOPPmB52\nAEzMsdsbsX37AVy+PATvv++C273z+h5lT7wRSUmvwGodh1tu+UZYzkDq9f/P//wFbnetbI/S8Xix\n2dajvv75wX8ohjEg7ACYuOIfEpIb5cCwzNixJRg3bhhGjBiNYcN6fQ5BHtoZNqwX99xzC5qbz6Gz\n8zN88okF3d07oN7j57APIzZabGdUJ4IxYuA/YUz+SCknjzWiq2ssurr6RgfvvPMKRo9+FZ9+Ovy6\nkfe+f/Dga+jt3QF/hyKfmQwAXiPPM3wZJjzYATCakU8YO3LkY7jd0h7l4yV3CN7ee3f3Lpw5owzl\nOK4bf+U51PT99aiqWsFGn2HCgB0AExGkCWPeRK00GlD22IONDpSPoXxbfg7JyK9HSspZfPObt3KP\nn2E0wA6AiSjy0UBHx2f45JMnZaEduTEPZuSV28pe/9zrvf4SNvwMoxFOAjNRRa4WunSpA+fPj0BX\n138gdLK4EQkJr8nCQH1KotTU4SzvZBgZrAJiDIPkEDo6PpWpewClkb/77nFoaTmPy5eHXE/ustFn\nGDXYATCGRD46YCPPMOHBDoBhGEZQtNjOGyLcFoZhGMYgsANgGIYRFHYADMMwgsIOgGEYRlDYATAM\nwwgKOwCGYRhBYQfAMAwjKOwAGIZhBIUdAMMwjKCwA2AYhhEUdgAMwzCCwg6AYRhGUNgBMAzDCAo7\nAIZhGEFhB8AwDCMo7AAYhmEEhR0AwzCMoITtAH7/+98jKysLQ4YMwbFjx4Iel5GRgenTpyM3Nxff\n/OY3w70cwzAME2HCdgDTpk3Dn/70J8ydG3oNV4vFAqfTiffeew9HjhwJ93JC4XQ6490E3cD3og++\nF33wvYgMYTuAzMxMTJ48eUDH8lq/g4Mf7j74XvTB96IPvheRIeo5AIvFgm9/+9uYOXMmfvWrX0X7\ncgzDMMwASQi1s6CgAF1dXQHvV1ZWYsGCBQO6QFNTE8aNG4fPPvsMBQUFyMzMxJw5c8JrLcMwDBMx\nLKQxPnPffffhpz/9KWbMmNHvsc899xy+8Y1v4Omnnw5siMWipRkMwzDCEq4ZDzkC0Hrxf/zjH7h2\n7RqGDx+OL7/8Eg6HAxs2bBjUORiGYZjoEHYO4E9/+hPS09PR0tKC4uJiPPDAAwCAc+fOobi4GADQ\n1dWFOXPm4I477sCsWbMwf/58FBYWRqblDMMwjCY0h4AYhmEYYxL3mcD19fXIzMzEpEmT8OKLL8a7\nOTGlvb0d9913H7KyspCdnY3q6moAwBdffIGCggJMnjwZhYWF8Hg8cW5p7Lh27Rpyc3N9IgNR74XH\n48H3vvc93H777Zg6dSoOHz4s7L144YUXkJWVhWnTpuHRRx9FT0+PMPfi8ccfx5gxYzBt2jTfe6E+\n+wsvvIBJkyYhMzMTDoej3/PH1QFcu3YNP/rRj1BfX4+TJ0/i9ddfx6lTp+LZpJgydOhQ/OxnP8OH\nH36IlpYWvPLKKzh16hQ2b96MgoICnD59GvPmzcPmzZvj3dSYUVVVhalTp/pEAaLei9WrV6OoqAin\nTp3C+++/j8zMTCHvRVtbG371q1/h2LFj+OCDD3Dt2jXU1tYKcy+WL1+O+vp6v/eCffaTJ09i165d\nOHnyJOrr67FixQp89dVXoS9AceTQoUNks9l82y+88AK98MILcWxRfPnnf/5nOnDgAE2ZMoW6urqI\niOj8+fM0ZcqUOLcsNrS3t9O8efPo4MGDNH/+fCIiIe+Fx+OhCRMmBLwv4r24ePEiTZ48mb744gu6\nevUqzZ8/nxwOh1D3wuVyUXZ2tm872GevrKykzZs3+46z2WzU3Nwc8txxHQF0dnYiPT3dt52WlobO\nzs44tih+tLW14b333sOsWbNw4cIFjBkzBgAwZswYXLhwIc6tiw0/+clPsHXrVtxwQ99jKeK9cLlc\nuPnmm7F8+XLMmDEDP/zhD/Hll18KeS9uuukmPP3007j11ltxyy23IDk5GQUFBULeC4lgn/3cuXNI\nS0vzHTcQexpXB8Dafy9///vf8d3vfhdVVVUYPny43z6LxSLEfaqrq8Po0aORm5sbVBIsyr3o7e3F\nsWPHsGLFChw7dgxf//rXA0IcotyL1tZWvPTSS2hra8O5c+fw97//Hb/73e/8jhHlXqjR32fv777E\n1QGkpqaivb3dt93e3u7nwUTg6tWr+O53v4vHHnsMDz74IACvV5dmYJ8/fx6jR4+OZxNjwqFDh7Bn\nzx5MmDABjzzyCA4ePIjHHntMyHuRlpaGtLQ03HXXXQCA733vezh27BjGjh0r3L149913ce+992LU\nqFFISEjAQw89hObmZiHvhUSw34TSnnZ0dCA1NTXkueLqAGbOnImPP/4YbW1tuHLlCnbt2oWFCxfG\ns0kxhYhQUlKCqVOn4sc//rHv/YULF2Lnzp0AgJ07d/ocg5mprKxEe3s7XC4Xamtr8a1vfQu//e1v\nhbwXY8eORXp6Ok6fPg0AePvtt5GVlYUFCxYIdy8yMzPR0tKC7u5uEBHefvttTJ06Vch7IRHsN7Fw\n4ULU1tbiypUrcLlc+Pjjj/svwR/phMVg2bt3L02ePJmsVitVVlbGuzkx5Z133iGLxUI5OTl0xx13\n0B133EH79u2jixcv0rx582jSpElUUFBAbrc73k2NKU6nkxYsWEBEJOy9+N///V+aOXMmTZ8+nb7z\nne+Qx+MR9l68+OKLNHXqVMrOzqalS5fSlStXhLkXixcvpnHjxtHQoUMpLS2Nfv3rX4f87Js2bSKr\n1UpTpkyh+vr6fs/PE8EYhmEEJe4TwRiGYZj4wA6AYRhGUNgBMAzDCAo7AIZhGEFhB8AwDCMo7AAY\nhmEEhR0AwzCMoLADYBiGEZT/D72T+rREex4RAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xb155298c>"
       ]
      }
     ],
     "prompt_number": 76
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print len(rk2ylist), len(rk4ylist)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "101 101\n"
       ]
      }
     ],
     "prompt_number": 83
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "temp = np.asarray(rk2ylist) - np.asarray(rk4ylist)\n",
      "plt.plot(temp)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 84,
       "text": [
        "[<matplotlib.lines.Line2D at 0xaeeab0ac>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtYVOXaBvB7FMpQ81SiAm1MQEARIZQy0TEdSRQytTyU\nmppSZqnZbrvbuy8oj9nJsm1maaaFmGVaKqnZmCeiBN0mbDGTRBBSkfJUKK7vj6fIAyeZGd41s+7f\ndc0lg8u1HlaxnnlPz2vSNE0DEREZTh3VARARkRpMAEREBsUEQERkUEwAREQGxQRARGRQTABERAZl\ncwIYPXo0PD09ERISUu7fW61WNGrUCGFhYQgLC8O0adNsvSQREdmBm60nGDVqFB5//HGMGDGiwmO6\nd++ONWvW2HopIiKyI5tbAFFRUWjSpEmlx3CtGRGR/jh8DMBkMmHHjh0IDQ1FTEwMMjMzHX1JIiKq\nBpu7gKoSHh6O3NxceHh4YP369ejfvz+ys7MdfVkiIqqCwxNAw4YNy77u06cPxo8fj6KiIjRt2vSy\n40wmk6NDISJySTXtZnd4F1BhYWFZcGlpadA07aqH/580TeNL0/Dcc88pj0EvL94L3gvei8pftrC5\nBTB06FBs2bIFx48fh4+PDxITE3H+/HkAQHx8PFauXIn58+fDzc0NHh4eWL58ua2XJCIiO7A5ASQl\nJVX694899hgee+wxWy9DRER2xpXAOmQ2m1WHoBu8F3/hvfgL74V9mDRbO5HsxGQy2dyfRURkNLY8\nO9kCICIyKCYAIiKDYgIgIjIoJgAiIoNiAiAiMigmACIig2ICICIyKCYAIiKDYgIgIjIoJgAiIoNi\nAiAiMigmACIig2ICICIyKCYAIiKDYgIgIjIoJgAiIoNiAiAiMigmACIig2ICICIyKCYAIgfjVtek\nV26qAyByBZoG7NsH7NgBbN8O7NoFnDwJnD4tr8aNAT8/oE0bIDQUuPdeICBAddRkdDa3AEaPHg1P\nT0+EhIRUeMwTTzwBf39/hIaGIiMjw9ZLEunGzz8Dc+bIwzwuTh7+XboAS5cC33wD/PQT8PvvQFYW\n8OqrQJ8+wOHDQPfukgimTQNOnFD9U5BhaTb6+uuvtfT0dK19+/bl/v3atWu1Pn36aJqmaampqVpk\nZGS5x9khFKJa89NPmjZihKY1aqRpDz2kaTt2aNrFi9X/9xcuaNrXX2vamDGa1qyZpiUmatqvvzou\nXnJdtjw7bW4BREVFoUmTJhX+/Zo1azBy5EgAQGRkJIqLi1FYWGjrZYmU+PVX4JlngLAwwNcXyMkB\nFi8G7rgDMJmqf566dYGoKOCdd4DUVCA7W7qI3nuPYwZUexw+CJyXlwcfH5+y997e3jhy5IijL0tk\nd599BgQGAvn5wJ49QGKi9O3bys8PWLYM+OIL4OWXgaFDgeJi289LVJVaGQTWrvhIY6rgo1JCQkLZ\n12azGWaz2YFREVXP2bPAU08B69cDK1YAXbs65jodOwJpacDTT8vXH3wA3HmnY65FzstqtcJqtdrl\nXA5PAF5eXsjNzS17f+TIEXh5eZV77KUJgEgPvv8euP9+IDwc2L0baNTIsde74QbgjTeA3r1lptB/\n/gMMGuTYa5JzufLDcWJiYo3P5fAuoLi4OLz//vsAgNTUVDRu3Bienp6OviyRzdatA+66C5g6Vbpo\nHP3wv1RsLLBhAzBxIvDWW7V3XTIWm1sAQ4cOxZYtW3D8+HH4+PggMTER58+fBwDEx8cjJiYG69at\ng5+fH+rXr4/FixfbHDSRo73xBjBjBvDppzKtU4WOHYGvv5bWwLFjwL//fW0DzURVMWlXdtArYjKZ\nrhorIKptFy8CkycDmzYBn38OtG6tOiKgoACwWIAHHpDWCNGlbHl2ciUw0R9KS4H4eFm0tX27fWb4\n2EOLFjJDqEsX+fqhh1RHRK6CCYAI8vAfNUpW6X7xBdCggeqILteqFZCSApjNwM03A337qo6IXAG7\ngMjwLlwAhg8Hjh8HVq8GPDxUR1Sx1FQZIF67FujcWXU0pAe2PDtZDZQMTdOAceOkHs9nn+n74Q8A\nt98uq4cHDZKBYSJbMAGQYWmaLLrKygJWrQLq1VMdUfXcc48MCA8bJl1XRDXFBECG9eKLsrp37Vqg\nfn3V0VybF16Qh78Na4CIOAhMxrRokSyw2rYNaNpUdTTXzs0NSEoCbrtNuoViYlRHRM6Ig8BkOJs3\nA0OGAFu3Am3bqo7GNtu2yXjAnj0AF9gbky3PTiYAMpT9+4Fu3eTT8113qY7GPp55RsYxPvmEK4WN\niLOAiKrhxAmgXz8p8eAqD38AeO454OBBqR5KdC3YAiBDOH9eyil07iyDv64mIwOIjpY/Kyi2Sy6K\nXUBEVZg0SXbd+uwz2Y3LFSUmykKxdevYFWQk7AIiqsSyZVLY7YMPXPfhD8hYQEEBu4Ko+tgCIJeW\nkSHllL/6CmjfXnU0jpeaCgwYAGRm6qeYHTkWu4CIylFUJPPkZ8+WXb2MYuxYKWkxd67qSKg2MAEQ\nXeHiRSmZ4O8PvPKK6mhq1/HjQLt2UtW0Y0fV0ZCjcQyA6AovvSQPwtmzVUdS+266CZg2DRg/XhIh\nUUWYAMjlbNsmn/qTkwF3d9XRqDFmjDz8lyxRHQnpGbuAyKUcOwaEhwMLFrA+TloacO+9Mv3V2Yrd\nUfWxC4gI8ol35EgplWz0hz8gi96ioow3BkLVxxYAuYxXX5Vun61bjdv1c6UffwQ6dZJpoSwW55o4\nC4gMb9cuoE8f4JtvgNatVUejL08+CZw7B8yfrzoScgQmADK0U6ek33/aNGDwYNXR6E9RkZS93roV\nCAxUHQ3Zm9IxgJSUFAQGBsLf3x+zy5lzZ7Va0ahRI4SFhSEsLAzTpk2z9ZJEl5kwQUo88+FfvqZN\ngX/8Q15El7JpR7DS0lJMmDABmzZtgpeXFzp16oS4uDgEBQVddlz37t2xZs0amwIlKk9SkpQ/SE9X\nHYm+TZggK4O//VbGBIgAG1sAaWlp8PPzg6+vL9zd3TFkyBCsXr36quPYtUOO8NNPwMSJwIcfcppj\nVerVk2Jx//d/qiMhPbEpAeTl5cHHx6fsvbe3N/Ly8i47xmQyYceOHQgNDUVMTAwyMzNtuSQRANkQ\nffhw4KmnpN4PVW30aJkNtGOH6khIL2zqAjJVo+h4eHg4cnNz4eHhgfXr16N///7Izs4u99iEhISy\nr81mM8xmsy3hkQubNUs2Rn/qKdWROI/rrweefVZ2ENu4UXU0VFNWqxVWq9Uu57JpFlBqaioSEhKQ\nkpICAJg5cybq1KmDf1Qy2tS6dWvs2rULTZs2vTwQzgKiakpLA2JjZeqnt7fqaJzL+fMyI+i992Tg\nnJyfsllAEREROHDgAHJyclBSUoLk5GTExcVddkxhYWFZcGlpadA07aqHP1F1nTkDPPggMG8eH/41\n4e4u4wDPPac6EtIDm7qA3NzcMG/ePERHR6O0tBRjxoxBUFAQFixYAACIj4/HypUrMX/+fLi5ucHD\nwwPLly+3S+BkTFOmAHfcAdx3n+pInNeDDwLTpwNWK8BeVmPjQjByGp99BjzxBLBnD3DjjaqjcW7v\nvgusWCF7BpBz40pgcnmFhbK5yYoVUuCMbFNSArRpA3z6KWdROTtWAyWXpmlS337UKD787eW662QG\n1cyZqiMhldgCIN176y1g4UJg5055cJF9nDkD3HqrjAVcsXifnAi7gMhl7d8PdO3KQmaOMm0a8MMP\nMi2UnBMTALmk8+eBLl2k62f8eNXRuKaTJwE/P6ml9Le/qY6GaoJjAOSSnn8euPlm4NFHVUfiupo0\nAR5+GHj5ZdWRkApsAZAubdsmc/0zMoAWLVRH49ry84H27YGDByUhkHNhC4Bcyi+/SKG3t9/mw782\ntGolpTX+WL9JBsIWAOnOiBGAh4fM/qHasWcPEBMDHDrEmVbOhi0AchnLl8u+vuyTrl2hoUBwsNx/\nMg62AEg3Dh8GIiKA9eu5OlWFlBTZNnL3bqAald5JJ9gCIKdXWipFyp58kg9/VaKj5b/Dl1+qjoRq\nCxMA6cL06VKq+OmnVUdiXCaTJGB2vxkHu4BIuW3bgEGDZDFSq1aqozG233+XBWFWK1deOwt2AZHT\nOnlSun4WLuTDXw+uvx4YN0423CHXxxYAKaNpwP33A56efODoyZ8Lww4dAho1Uh0NVYUtAHJKb70l\nhcheekl1JHSpVq2Au+8GFi9WHQk5GlsApMTu3YDFAmzfDgQEqI6GrpSaCjzwAJCdDdStqzoaqgxb\nAORUTp2Srp+5c/nw16vISKBpU1mTQa6LLQCqVZomg74eHjLwS/q1dKm8NmxQHQlVhvsBkNOYP19e\nqamSBEi//pwSunmzlIkgfWICIKeQlgb06yf9/v7+qqOh6nj2WaC4GHjjDdWRUEWYAEj3jh+XEg9z\n5wL9+6uOhqorN1cKxR0+DDRooDoaKo/SQeCUlBQEBgbC398fs2fPLveYJ554Av7+/ggNDUVGRoat\nlyQnU1oqM0qGDOHD39n4+ADduwMffKA6EnIEmxJAaWkpJkyYgJSUFGRmZiIpKQlZWVmXHbNu3Tr8\n8MMPOHDgAN5++208yv39DOfZZ4GSEqn3Q85n/HjgP/+RAXxyLTYlgLS0NPj5+cHX1xfu7u4YMmQI\nVq9efdkxa9aswciRIwEAkZGRKC4uRmFhoS2XJSeSnAwkJQErVgBubqqjoZro2RP47Tdgxw7VkZC9\n2ZQA8vLy4OPjU/be29sbeXl5VR5z5MiRCs/5668ySLhvnyxJ/+03WyIkldLTgQkTgE8/lc3dyTnV\nqQM8+qi0Asi12PSZzFTNXSOuHKCo6N8lJCRg505JAM2amVFSYkZRkWxU3bq1vNq3B8LD5dW8uS3R\nkyP9/DNw770y5TM0VHU0ZKuRI4HERPnvyt87taxWK6xWq13OZVMC8PLyQm5ubtn73NxceHt7V3rM\nkSNH4OXlVe75EhIS8MILQKdOwLRp8r2LF4GjR6Uw1Y8/Anv3AnPmyKfLRo2AHj3kddddwBWXJkV+\n+00e/iNHSplncn5NmgADBwLvvAM884zqaIzNbDbDbDaXvU9MTKzxuWzqAoqIiMCBAweQk5ODkpIS\nJCcnIy4u7rJj4uLi8P777wMAUlNT0bhxY3h6elZ4znPngBtuuCTAOoCXF9C1q2wWPmeO7FhUVCTL\n1Dt1Aj77TD5lhocDL7wgSYIDVmpcvCgP/ltuARISVEdD9vToo7J6u7RUdSRkLza1ANzc3DBv3jxE\nR0ejtLQUY8aMQVBQEBYsWAAAiI+PR0xMDNatWwc/Pz/Ur18fi6soMXj2bPX6i00mIChIXuPHAxcu\nSNfRp58CsbGyynTECJl+eMkQBDnYM8/I2M3GjZK8yXXcdhvQrJmUhujTR3U0ZA+6WwgWHw+EhQGP\nPFLzc2mazFh4/31g5Ur5H/eRR4C4OM5EcaQFC2Q7wZ075UFBrmfhQmDtWvmgRfrgUtVAz52zvUaM\nyQTceac8kPLypEvilVcAX1/g+edlIIvs65NPpMtn3To+/F3Z0KHA11/L7xU5P90lgLNnLx8DsFW9\netINtG2bPJyOHAHatpUWwf799ruOkW3aJPdz7VrAz091NORIDRrIiu5331UdCdmD7hLAlYPA9tSh\nA/D22/Lg9/QEoqJkZkN6umOuZwSpqfKpcOVKGYQn1xcfL11BFy6ojoRspbsEcPas48sEN28uc5oP\nHZIkEBsrVSpTUx17XVezZw9wzz3AkiVAt26qo6HaEhoqU665WYzz010CcGQL4Er16wOTJgEHDwJ9\n+wKDB8uf331XO9d3ZunpQHS0bOYeE6M6Gqpt8fEyxkbOTZcJoLY3CqlXT+Y4Z2dLAujfXz7Z/ve/\ntRuHs/j2W5kGOH8+cN99qqMhFe6/X2Z7HT6sOhKyhe4SgL0Hga/F9dfLmoIDBwCzGejdWwaQf/hB\nTTx6lJoqSfKdd2S1LxmTh4eM/SxapDoSsoXuEkBtdgFV5IYbgMmTJREEBQG33y4thPx8tXGp9vnn\nspZiyRIZNyFjGztWZgNxZbDz0mUC0MtesQ0bAv/+t8waatgQCAkBpk4FTp5UHVntW7BAfuE//5yr\nQEmEhgKtWgEpKaojoZrSXQJQ2QVUkWbNgBdflFkvRUVAQAAwcyZw5ozqyBzv4kUp7/DSS8DWrUDn\nzqojIj0ZN06mVpNz0lUpiNJSDW5uMr9Yz3VksrNll6utW4F//Qt4+GEZP3A1J08Cw4fLpuCrVrGm\nP13t9Gkp/Ld3rxRtpNrnMqUgfvtNHqR6fvgD0gJITpbukLVrZWXxu+8C58+rjsx+9uwBIiJkZe9X\nX/HhT+Vr0EBmglVR45F0SlePWj0MAF+L8HApL/Hhh/IKCpIBUmdOBJoGvPUW0KuX7Mnw2muAu7vq\nqEjPxo2TWWEXL6qOhK6V7hKAXgaAr0WXLrJHwcKFwHvvSYtg4ULZCN2Z5ObK4q5Fi6Tg19ChqiMi\nZ3DbbUDTplICnJyLrhKAHgeAr0WPHtJd8mcZ6jZtZPD0119VR1a50lIZyAsPl/UPO3ZIa4aousaO\nlQ895Fx0lQCctQVwpa5dgS++kJrpu3bJXsZPPw3k5KiO7GpffSV9/e+/L62YZ57hngl07YYNk6qw\nLLXuXHSVAJy9BXCl224DkpIkCVy4IA/afv1k4Fj14pnvvpNyF6NHy0N/61aplkpUE40aycrwJUtU\nR0LXQlcJwNkGgavL11c2pDl8WMpPJybK1LmnngJ27669/Ys1TSo43nUXMGCAdPdkZcksDpOpdmIg\n1zV2rAwG62NiOVWH7hKAK3QBVcTDAxg1CkhLkwGz66+XwnPBwdJFtGWLY2YQ7d0L/POfwK23yp+j\nR0sF1MmTpRAekT3ccQdQt660Jsk56Goh2IoVGpKTZQDVKC5elOqaa9fK69Ah+UW64w6pQRQWdm1b\nLGqazObZulVm8litkliHDpVXaCg/7ZPjvPqqlApfulR1JMZhy0IwXSWAJUs0bNokA5JGVVAgs3B2\n7pTX3r3yqcrPT7qSGjcGbrxRahNduCDjJmfPAoWF8qn+4EFZnHPnnUD37rLhTWio/hfXkWs4flz+\nXz10CGjSRHU0xuAyCWD+fA0ZGdxo4lKaBpw4IZVJc3JkSumpU/Knm5t0K3l4ADfdJL94bdrIgByR\nKkOHytqYxx9XHYkx2JIAdDXhz1UHgW1hMsnD/aabpFuISO/GjpWd9iZMYHej3tW4Y6CoqAgWiwUB\nAQHo3bs3iouLyz3O19cXHTp0QFhYGDpXUUrS1QeBiYzAbJZKud9+qzoSqkqNE8CsWbNgsViQnZ2N\nnj17YtasWeUeZzKZYLVakZGRgbS0tErP6WrrAIiMqE4dqZDLlcH6V+MEsGbNGowcORIAMHLkSHz6\n6acVHlvd/im2AIhcw0MPyWy+U6dUR0KVqXECKCwshKenJwDA09MThYWF5R5nMpnQq1cvREREYGEV\nHwnYAiByDS1byiy05GTVkVBlKh0EtlgsKCgouOr706dPv+y9yWSCqYLRnu3bt6Nly5Y4duwYLBYL\nAgMDERUVVe6xO3cmID9f6omYzWaYzeZq/hhEpDdjxwLPPy/dQWQ/VqsVVqvVLueq8TTQwMBAWK1W\ntGjRAkePHkWPHj3wv//9r9J/k5iYiAYNGmDKlClXB2Iy4f77NQwYAAweXJOIiEhPSktl7cratawz\n5UhKdgSLi4vDkj8qPy1ZsgT9+/e/6pizZ8/i1B+dgGfOnMGGDRsQEhJS4TnZBUTkOurWldIn77yj\nOhKqSI0TwNSpU7Fx40YEBARg8+bNmDp1KgAgPz8fffv2BQAUFBQgKioKHTt2RGRkJPr164fevXtX\neE4OAhO5ljFjgA8+kN9t0h9drQS+4w4Nc+ZIGQMicg19+sh+AcOHq47ENbnMpvBcCUzkerhbmH7p\nLgGwC4jItcTGSi2rrCzVkdCVdJUAOAhM5Hrc3WUwmK0A/dHVGMBNN2nIzARuvll1NERkTz/+CERG\nyl4V3ITIvlxmDIAtACLXdOutQMeOwKpVqiOhS+kqAXAQmMh1jRsHvP226ijoUrrqArruOg2//646\nEiJyhJISwMdHtisNCFAdjetwmS4gfvoncl3XXSdVQtkK0A9dtQBattSQn686EiJylIMHZWe7w4c5\nGGwvbAEQkVNo0wYICwM+/lh1JK7B1o/vTABEVKvi44G33lIdhfP77TegfXvbzqGrBMBVwESuLzZW\nuoL27VMdiXP74gvb10zpKgGwBUDk+tzdZZOYBQtUR+LckpNt3ztFV4PAd9+tYf161ZEQkaMdPixj\nAYcPA/Xrq47G+Zw9C7RqBWRnA56eHAQmIidyyy1S9p17BtfM+vVARATQvLlt52ECICIlHnkEmD9f\ndRTOyR7dP4DOEgAHgYmMIzoaOHECSEtTHYlzOX1aBoAHDLD9XLpKAGwBEBlH3brAo48Cb76pOhLn\nsnatLKZr1sz2c7nZfgr7YQuAyFhGj5bFYcePAzfdpDoa52Cv7h+ALQAiUqhZM6B/f+Ddd1VH4hx+\n/RX48ku5Z/bABEBESj32mAwGl5aqjkT/Vq8GoqKAJk3scz5dJQB2AREZT6dOMp2Ra4CqtmwZ8OCD\n9jufrhIAWwBExvTYYxwMrsrRozJj6p577HfOGieAjz76CO3atUPdunWRnp5e4XEpKSkIDAyEv78/\nZs+eXek52QIgMqbBg4Fdu2RlK5Xvww+Be++17wflGieAkJAQrFq1Ct26davwmNLSUkyYMAEpKSnI\nzMxEUlISsrKyKjyeLQAiY6pXDxg7FnjjDdWR6NfSpfbt/gFsSACBgYEIqGJft7S0NPj5+cHX1xfu\n7u4YMmQIVq9eXeHxTABExjV+PPDBB0BxsepI9GfvXlk0Zzbb97wOHQPIy8uDj49P2Xtvb2/k5eVV\neDy7gIiMy8sLuPtuYNEi1ZHoz7JlwAMPAHXs/MSudCGYxWJBQUHBVd+fMWMGYmNjqzy5yWS6pmCW\nLk3Al1/K12azGWZ7pzsi0rWJE4EhQ+TPunVVR6MPpaXSMvriC3lvtVphtVrtcu5KE8DGjRttOrmX\nlxdyc3PL3ufm5sLb27vC4ydNSkBwsE2XJCInFhkJtGgBrFkjA54EWK0yTbZdO3l/5YfjxMTEGp/b\nLg2KimpRR0RE4MCBA8jJyUFJSQmSk5MRFxdX4Xk4BkBEEycCc+eqjkI/HDH4+6caJ4BVq1bBx8cH\nqamp6Nu3L/r06QMAyM/PR9++fQEAbm5umDdvHqKjoxEcHIzBgwcjKCiownMyARDRwIHADz8Au3er\njkS9X36R1b+OSgC62hHsl1803Hij6kiISLXZs4Hvv5dPv0Y2fz6weTPw0UcVH2My1XxHMF0lgJIS\nDe7uqiMhItWKi6VKaEaG7B5mVOHhwKxZQO/eFR9jSwLQVSkIPvyJCAAaNwZGjQJee011JOrs2gUU\nFQG9ejnuGrpqAegkFCLSgSNHgA4dgIMH7Vf90pk88oisjXj22cqPc5kuIJ2EQkQ6MXIkEBgI/POf\nqiOpXadPAz4+sgK4kpnzAFyoC4iI6FJPPQW8/jrw22+qI6ldK1YAXbtW/fC3FRMAEelWSAgQFial\nEIxk4UIpjudo7AIiIl3bsgV4+GEgKwtw09Uu5o6xezfQrx+Qk1O9n5ddQETksrp1k/IQycmqI6kd\nc+fKBjm1kezYAiAi3du4UUpE7N3r2kXifv4ZaNtWVkI3a1a9f8MWABG5tF69gEaNgI8/Vh2JYy1Y\nAAwaVP2Hv63YAiAip7BuHfCPfwB79ti/Lr4elJQAf/sbsGGDDH5XF1sAROTy+vQBrr9eiqO5ohUr\ngODga3v424oJgIicgskkq2Kffx5wtc4CTZPB30mTave6TABE5DTi4iQRuNpYwM6dwMmTwB+V9GsN\nEwAROQ2TCZg5E/jXv4ALF1RHYz8zZgCTJ9f+2AYTABE5ld69pUja4sWqI7GPjAx5jRlT+9fmLCAi\ncjppacCAAUB2NuDhoToa2wwcKHV/Jk+u2b9nNVAiMpxBg4BOnWRqqLP6/ntZ4/DjjzVPZEwARGQ4\n+/fLJ+fsbOfdL2DoUKBjR9uSGBMAERlSfLx8cn71VdWRXLv9+4GoKNnwpmHDmp+HCYCIDOnYMVk8\nZbUC7dqpjubajBgB+PtXveNXVZgAiMiw5s0DPvkE+PJLmSbqDDIyZGXz/v1S48gWLAVBRIb1yCPA\niRPARx+pjqR6NA2YMgV47jnbH/62qnEC+Oijj9CuXTvUrVsX6enpFR7n6+uLDh06ICwsDJ07d67p\n5YiIyuXmBrzxhmwfeeaM6miq9vnnQEFB7ez4VZUaJ4CQkBCsWrUK3bp1q/Q4k8kEq9WKjIwMpKWl\n1fRyREQV6tZNBlSnTVMdSeXOnwf+/nfgpZf0sbtZjUMIDAys9rHs2yciR5szR6ZU3ncfEB6uOpry\nLVgA3HKL9P/rgcPHAEwmE3r16oWIiAgsXLjQ0ZcjIoNq1Qp45RWZXfP776qjudqJE8ALL8inf70M\nVlfaArBYLCgoKLjq+zNmzEBsbGy1LrB9+3a0bNkSx44dg8ViQWBgIKKioso9NiEhoexrs9kMs9lc\nrWsQEQHAAw8AK1cCiYlSYE1PJk6UhV8dOth2HqvVCqvVapeYbJ4G2qNHD7z88ssIr0abKzExEQ0a\nNMCUKVOuDoTTQInIDgoLgdBQ2TgmMlJ1NGLNGuDJJ2U3s/r17Xtu5dNAK7r42bNncerUKQDAmTNn\nsGHDBoTU5nY3RGQ4np7A668DI0cCZ8+qjgYoKgIefRRYtMj+D39b1TgBrFq1Cj4+PkhNTUXfvn3R\n549Rjfz8fPT9Y1eDgoICREVFoWPHjoiMjES/fv3Qu3dv+0RORFSB+++XT/9jx6rfPWzyZKn4WcWE\nSSW4EpiIXNK5c8CddwLDh9e81LKtVq+Wa+/d67hP/ywFQURUjpwc4PbbgaQkoEeP2r12VhbQvbv0\n/99+u+Ouo3wMgIhIj3x9gWXLgGHDgMOHa++6J0/K/sUvvujYh7+tmACIyKX16iX19i0W4OhRx1/v\nwgVg8GCVa2uVAAAHpElEQVSgXz/goYccfz1b6GAxMhGRY02aBJw+DfTsCXz1lcwUcgRNk+megKxM\n1jsmACIyhH//W2rx9OoFbN4M3Hyzfc9fWgpMmAB8+y2wcaM+av1UxQlCJCKyj4QESQJms+wh0Lat\nfc57/ryUoDh6VJLLjTfa57yOxjEAIjIMkwmYPl3KMnTtap89BIqLgf79pRT1+vXO8/AHmACIyGBM\nJmDcOCAlRQaHJ06s+T4Cq1bJVpStWwMffwzccIN9Y3U0JgAiMqTbbgN27ZJ9hdu0kUHb6iaCH3+U\n1b1Tp8oag3nzAHd3x8brCEwARGRYTZoAH34IbNokg7dt2sgsnpUrgbw8OUbTZFXxTz8Bb74pq4sj\nI+WT/549+izxUF1cCUxE9Id9+6R8w86d8rpwQQrK1akjyaJnT1lUZrHo5xM/S0EQEdmZpskmLg0a\nAPXqqY6mYkwAREQGxVpARER0zZgAiIgMigmAiMigmACIiAyKCYCIyKCYAIiIDIoJgIjIoJgAiIgM\nigmAiMigapwA/v73vyMoKAihoaEYMGAAfvnll3KPS0lJQWBgIPz9/TF79uwaB0pERPZV4wTQu3dv\n7Nu3D3v27EFAQABmzpx51TGlpaWYMGECUlJSkJmZiaSkJGRlZdkUsBFYrVbVIegG78VfeC/+wnth\nHzVOABaLBXXqyD+PjIzEkSNHrjomLS0Nfn5+8PX1hbu7O4YMGYLVq1fXPFqD4P/cf+G9+AvvxV94\nL+zDLmMAixYtQkxMzFXfz8vLg4+PT9l7b29v5P1ZZJuIiJSqdFN4i8WCgoKCq74/Y8YMxMbGAgCm\nT5+O6667DsOGDbvqOJPJZKcwiYjI7jQbLF68WOvSpYt27ty5cv9+586dWnR0dNn7GTNmaLNmzSr3\nWAB88cUXX3zV4FVTlbYAKpOSkoI5c+Zgy5YtqFfBbgkRERE4cOAAcnJy0KpVKyQnJyMpKancYzXu\nBUBEVKtqPAbw+OOP4/Tp07BYLAgLC8P48eMBAPn5+ejbty8AwM3NDfPmzUN0dDSCg4MxePBgBAUF\n2SdyIiKyiW52BCMiotqlfCWwkReK5ebmokePHmjXrh3at2+P119/HQBQVFQEi8WCgIAA9O7dG8XF\nxYojrT2lpaUICwsrm2Rg1HtRXFyMQYMGISgoCMHBwfjmm28Mey9mzpyJdu3aISQkBMOGDcPvv/9u\nmHsxevRoeHp6IiQkpOx7lf3sM2fOhL+/PwIDA7Fhw4Yqz680ARh9oZi7uzteffVV7Nu3D6mpqXjz\nzTeRlZWFWbNmwWKxIDs7Gz179sSsWbNUh1pr5s6di+Dg4LIZZEa9FxMnTkRMTAyysrLw3//+F4GB\ngYa8Fzk5OVi4cCHS09Oxd+9elJaWYvny5Ya5F6NGjUJKSspl36voZ8/MzERycjIyMzORkpKC8ePH\n4+LFi5VfoMbDx3awY8eOy2YJzZw5U5s5c6bCiNS65557tI0bN2pt27bVCgoKNE3TtKNHj2pt27ZV\nHFntyM3N1Xr27Klt3rxZ69evn6ZpmiHvRXFxsda6deurvm/Ee3HixAktICBAKyoq0s6fP6/169dP\n27Bhg6HuxaFDh7T27duXva/oZ79ylmV0dLS2c+fOSs+ttAXAhWJ/ycnJQUZGBiIjI1FYWAhPT08A\ngKenJwoLCxVHVzsmT56MOXPmlK0wB2DIe3Ho0CHcfPPNGDVqFMLDwzF27FicOXPGkPeiadOmmDJl\nCm655Ra0atUKjRs3hsViMeS9+FNFP3t+fj68vb3LjqvO81RpAuBCMXH69GkMHDgQc+fORcOGDS/7\nO5PJZIj79Pnnn6N58+YICwurcEqwUe7FhQsXkJ6ejvHjxyM9PR3169e/qovDKPfi4MGDeO2115CT\nk4P8/HycPn0ay5Ytu+wYo9yL8lT1s1d1X5QmAC8vL+Tm5pa9z83NvSyDGcH58+cxcOBADB8+HP37\n9wcgWf3PFdhHjx5F8+bNVYZYK3bs2IE1a9agdevWGDp0KDZv3ozhw4cb8l54e3vD29sbnTp1AgAM\nGjQI6enpaNGiheHuxXfffYcuXbqgWbNmcHNzw4ABA7Bz505D3os/VfQ7ceXz9MiRI/Dy8qr0XEoT\nwKULxUpKSpCcnIy4uDiVIdUqTdMwZswYBAcHY9KkSWXfj4uLw5IlSwAAS5YsKUsMrmzGjBnIzc3F\noUOHsHz5ctx1111YunSpIe9FixYt4OPjg+zsbADApk2b0K5dO8TGxhruXgQGBiI1NRXnzp2DpmnY\ntGkTgoODDXkv/lTR70RcXByWL1+OkpISHDp0CAcOHEDnzp0rP5m9Byyu1bp167SAgACtTZs22owZ\nM1SHU6u2bt2qmUwmLTQ0VOvYsaPWsWNHbf369dqJEye0nj17av7+/prFYtFOnjypOtRaZbVatdjY\nWE3TNMPei927d2sRERFahw4dtHvvvVcrLi427L2YPXu2FhwcrLVv314bMWKEVlJSYph7MWTIEK1l\ny5aau7u75u3trS1atKjSn3369OlamzZttLZt22opKSlVnp8LwYiIDEr5QjAiIlKDCYCIyKCYAIiI\nDIoJgIjIoJgAiIgMigmAiMigmACIiAyKCYCIyKD+H0sufbAOJ3l3AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xb1432c0c>"
       ]
      }
     ],
     "prompt_number": 84
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "anharmonic oscillator\n",
      "===\n",
      "$$-m\\frac{d^2x}{dt^2} -10\\frac{x}{|x|} +5\\cos(\\omega t) = 0$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}