guyer/newton.ipynb

## newton.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform Newton iterations to solve non-linear equations in FiPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import fipy as fp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mesh = fp.Grid1D(nx=100, dx=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## single equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "C = fp.CellVariable(mesh=mesh, name=\"C\", hasOld=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "fipy/viewers/matplotlibViewer/__init__.py:113: UserWarning: Matplotlib1DViewer efficiency is improved by setting the 'datamax' and 'datamin' keys\n",
      "  return Matplotlib1DViewer(vars=vars, title=title, axes=axes, **kwlimits)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARcAAAEKCAYAAAA4mxGRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEf1JREFUeJzt3X+s3Xd93/HnCzuh7chiEjQHx67CiqPZbHShaeas27jl\nR3vjdjYIqamrFppKA20JRVnHnBA0PLSJZT9IilJoBAG8allWQQGzJSQe5e4PlAaypiklNok33NmJ\n4kCqIBbE6sjv/XG+RCeHc+917tcfn+91ng/pKOf7+Xy+3/O+0Tmv+/l+v597nKpCkk61F826AEln\nJsNFUhOGi6QmDBdJTRgukpowXCQ1YbhIasJwURNJfiXJ/Um+m+SxJHcm+ZlZ16XTx3DRKZfknwI3\nAf8K+GvAJuB3gB2zrEunV1yhq1MpybnAUeDXq+rTs65Hs+PMRafa5cCPAJ+ZdSGaLcNFp9r5wLer\n6sSsC9FsGS461Z4EXpbE99YLnG8AnWr3Av8PePOsC9FsGS46parqO8C/AH4nyc4kP5bkrCRXJLlx\n1vXp9PFukZpI8ivAtcAW4LvA/cC/rqo/mmlhOm16h0uSeeBmYA3wsar6od9OST4EXAF8j9Etyge6\n9nXAx4BXAQX8hm8+6czQ67QoyRrgFmAe2ArsSrJlYsx24JVVtRl4O/CRse7fBu6sqi3Aq4EDfeqR\nNBx9r7lcBhyqqsNVdRy4A9g5MWYHsBegqu4D1iVZ3y22+vtV9fGu75nufF3SGaBvuFwIHBnbPtq1\nLTdmI/AK4FtJPpHkj5N8NMmP9axH0kD0DZeTvWCTKfutBV4DfLiqXgM8DVzXsx5JA7G25/6PMvqj\ntB/YxGhmstSYjV1bgKNV9dWu/VNMCZck3s6SZqSqJicGJ63vzOV+YHOSi5KcDVwJ7JsYsw94K0CS\nbcBTVXWsqh4HjiS5uBv3BuDr016kqlbV433ve9/MaziT67Xm0/Poq9fMpaqeSXINcDejW9G3VdWB\nJO/o+m+tqjuTbE9yiNGpz1Vjh3gn8J+6YPpfE32SVrG+p0VU1V3AXRNtt05sX7PIvg8CP923BknD\n4/L/Bubm5mZdwvOy2uoFa14NBr/8P0kNvUbpTJSE6nFBt/dp0awkK/6ZB8fw1Jlo1YYLnBkfyjMp\nJKVxXnOR1IThIqkJw0VSE4aLpCYMl4Zuv/12Lr30Us455xw2bNjA9u3b+fKXvzzrsqTTwnBp5IMf\n/CDXXnst733ve3niiSc4cuQIV199Nfv2Tf7plXRmWrWL6LoFPjOoaHnf+c532LhxI5/85Cd5y1ve\nsuTYIf8cemHru4jOmUsD9957L9///vd585v91zX0wrWqF9Et51SsT1vJpOLJJ5/kZS97GS96kdmt\nF64zOlxmdbZx/vnn8+1vf5sTJ04YMHrB8p3fwOWXX86LX/xiPvMZ/y12vXAZLg2ce+65vP/97+fq\nq6/mc5/7HN/73vc4fvw4d911F7t37551edJp4d2ihm6//XZuuukmDhw4wDnnnMOll17KDTfcwLZt\n254dsxp+Dr0w9b1bZLjM2Jnyc+jM461oSYNkuEhqwnCR1IThIqkJw0VSE4aLpCZW9fJ/v9xaGq5V\nGy6uDZGGzdMiSU0YLpKa6B0uSeaTHEzySJKpf5WX5ENd/4NJLpnoW5PkgSSf71uLpOHoFS5J1gC3\nAPPAVmBXki0TY7YDr6yqzcDbgY9MHOZdwEOAF1GkM0jfmctlwKGqOlxVx4E7gJ0TY3YAewGq6j5g\nXZL1AEk2AtuBjwHe+pHOIH3D5ULgyNj20a7tZMfcBLwbONGzDkkD0zdcTvZUZnJWkiS/CDxRVQ9M\n6Ze0yvVd5/IosGlsexOjmclSYzZ2bW8BdnTXZH4E+KtJ/mNVvXXyRfbs2fPs87m5Oebm5nqWLWnS\nwsICCwsLp+x4vb4sKsla4BvA64HHgK8Au6rqwNiY7cA1VbU9yTbg5qraNnGc1wL/rKr+4ZTXmPpl\nUZLa6vtlUb1mLlX1TJJrgLuBNcBtVXUgyTu6/lur6s4k25McAp4GrlrscH1qkTQsq/ZrLiW15ddc\nShokw0VSE4aLpCYMF0lNGC6SmjBcJDVhuEhqwnCR1IThIqkJw0VSE4aLpCYMF0lNGC6SmjBcJDVh\nuEhqwnCR1IThIqkJw0VSE4aLpCYMF0lNGC6SmjBcJDVhuEhqwnCR1IThIqkJw0VSE4aLpCYMF0lN\nGC6SmugdLknmkxxM8kiS3YuM+VDX/2CSS7q2TUm+lOTrSf4syW/2rUXScPQKlyRrgFuAeWArsCvJ\nlokx24FXVtVm4O3AR7qu48C1VfUqYBtw9eS+klavvjOXy4BDVXW4qo4DdwA7J8bsAPYCVNV9wLok\n66vq8ar6k679/wIHgA0965E0EH3D5ULgyNj20a5tuTEbxwckuQi4BLivZz2SBmJtz/3rJMdlsf2S\nvAT4FPCubgbzQ/bs2fPs87m5Oebm5p5XkZKWt7CwwMLCwik7XqpONh+m7JxsA/ZU1Xy3fT1woqpu\nHBvzu8BCVd3RbR8EXltVx5KcBfxX4K6qunmR16g+NUpamSRU1eTE4KT1PS26H9ic5KIkZwNXAvsm\nxuwD3grPhtFTXbAEuA14aLFgkbR69TotqqpnklwD3A2sAW6rqgNJ3tH131pVdybZnuQQ8DRwVbf7\nzwC/Cvxpkge6tuur6gt9apI0DL1Oi04HT4uk2Zj1aZEkTWW4SGrCcJHUhOEiqQnDRVIThoukJgwX\nSU0YLpKaMFwkNWG4SGrCcJHUhOEiqQnDRVIThoukJgwXSU0YLpKaMFwkNWG4SGrCcJHUhOEiqQnD\nRVIThoukJgwXSU0YLpKaMFwkNWG4SGrCcJHUhOEiqYne4ZJkPsnBJI8k2b3ImA91/Q8mueT57Ctp\ndeoVLknWALcA88BWYFeSLRNjtgOvrKrNwNuBj5zsvpJWr74zl8uAQ1V1uKqOA3cAOyfG7AD2AlTV\nfcC6JBec5L6SVqm+4XIhcGRs+2jXdjJjNpzEvpJWqbU996+THJc+L5LsGdua6x6STq2F7nFq9A2X\nR4FNY9ubGM1AlhqzsRtz1knsC0DVnp5lSlreHOO/uJN/2etofU+L7gc2J7koydnAlcC+iTH7gLcC\nJNkGPFVVx05yX0mrVK+ZS1U9k+Qa4G5gDXBbVR1I8o6u/9aqujPJ9iSHgKeBq5bat089koYjVSd7\n2WQ2ktTQa5TOREmoqhVfL3WFrqQmDBdJTRgukpowXCQ1YbhIasJwkdSE4SKpCcNFUhOGi6QmDBdJ\nTRgukpowXCQ1YbhIasJwkdSE4SKpCcNFUhOGi6QmDBdJTRgukpowXCQ1YbhIasJwkdSE4SKpCcNF\nUhOGi6QmDBdJTRgukpowXCQ10StckpyXZH+Sh5Pck2TdIuPmkxxM8kiS3WPt/y7JgSQPJvmDJOf2\nqUfScPSduVwH7K+qi4EvdtvPkWQNcAswD2wFdiXZ0nXfA7yqqn4SeBi4vmc9kgaib7jsAPZ2z/cC\nb5oy5jLgUFUdrqrjwB3AToCq2l9VJ7px9wEbe9YjaSD6hsv6qjrWPT8GrJ8y5kLgyNj20a5t0m8A\nd/asR9JArF1uQJL9wAVTum4Y36iqSlJTxk1rm3yNG4C/rKrbp/Xv2bPn2edzc3PMzc0td0hJz9PC\nwgILCwun7HipWvazv/jOyUFgrqoeT/Jy4EtV9TcmxmwD9lTVfLd9PXCiqm7stn8d+EfA66vq+1Ne\no/rUKGllklBVWen+fU+L9gFv656/DfjslDH3A5uTXJTkbODKbj+SzAPvBnZOCxZJq1ffmct5wO8D\nPw4cBn6pqp5KsgH4aFX9QjfuCuBmYA1wW1V9oGt/BDgb+IvukPdW1T+ZeA1nLtIM9J259AqX08Fw\nkWZj1qdFkjSV4SKpCcNFUhOGi6QmDBdJTRgukpowXCQ1YbhIasJwkdSE4SKpCcNFUhOGi6QmDBdJ\nTRgukpowXCQ1YbhIasJwkdSE4SKpCcNFUhOGi6QmDBdJTRgukpowXCQ1YbhIasJwkdSE4SKpCcNF\nUhOGi6QmVhwuSc5Lsj/Jw0nuSbJukXHzSQ4meSTJ7in9v5XkRJLzVlqLpOHpM3O5DthfVRcDX+y2\nnyPJGuAWYB7YCuxKsmWsfxPwRuDPe9QhaYD6hMsOYG/3fC/wpiljLgMOVdXhqjoO3AHsHOv/IPDP\ne9QgaaD6hMv6qjrWPT8GrJ8y5kLgyNj20a6NJDuBo1X1pz1qkDRQa5fqTLIfuGBK1w3jG1VVSWrK\nuGltJPlR4D2MTomebV66VEmryZLhUlVvXKwvybEkF1TV40leDjwxZdijwKax7U2MZi8/AVwEPJgE\nYCPwP5NcVlU/dJw9e/Y8+3xubo65ubmlypa0AgsLCywsLJyy46Vq6uRi+R2Tfws8WVU3JrkOWFdV\n102MWQt8A3g98BjwFWBXVR2YGPdN4Keq6i+mvE6ttEZJK5eEqlrxGUWfay7/BnhjkoeB13XbJNmQ\n5L8BVNUzwDXA3cBDwH+ZDJaO6SGdYVY8czldnLlIszHLmYskLcpwkdSE4SKpCcNFUhOGi6QmDBdJ\nTRgukpowXCQ1YbhIasJwkdSE4SKpCcNFUhOGi6QmDBdJTRgukpowXCQ1YbhIasJwkdSE4SKpCcNF\nUhOGi6QmDBdJTRgukpowXCQ1YbhIasJwkdSE4SKpCcNFUhMrDpck5yXZn+ThJPckWbfIuPkkB5M8\nkmT3RN87kxxI8mdJblxpLZKGp8/M5Tpgf1VdDHyx236OJGuAW4B5YCuwK8mWru9ngR3Aq6vqbwL/\nvkctg7KwsDDrEp6X1VYvWPNq0CdcdgB7u+d7gTdNGXMZcKiqDlfVceAOYGfX94+BD3TtVNW3etQy\nKKvtTbTa6gVrXg36hMv6qjrWPT8GrJ8y5kLgyNj20a4NYDPwD5L8UZKFJJf2qEXSwKxdqjPJfuCC\nKV03jG9UVSWpKeOmtY2/9kuraluSnwZ+H/jry9QrabWoqhU9gIPABd3zlwMHp4zZBnxhbPt6YHf3\n/C7gtWN9h4DzpxyjfPjwMZvHSvOhqpaeuSxjH/A24Mbuv5+dMuZ+YHOSi4DHgCuBXV3fZ4HXAf8j\nycXA2VX15OQBqio9apQ0I+lmB89/x+Q8RqcyPw4cBn6pqp5KsgH4aFX9QjfuCuBmYA1wW1V9oGs/\nC/g48LeBvwR+q6oWev00kgZjxeEiSUsZ9ArdpRbgDUGSTUm+lOTr3ULA3+zaT2qB4SwlWZPkgSSf\n77YHXXOSdUk+1S26fCjJ3xlyzUmu7d4TX0tye5IXD63eJB9PcizJ18baFq0xyfXdZ/Fgkp9b7viD\nDZelFuANyHHg2qp6FaOL11d3NS67wHAA3gU8xOjCHQy/5t8G7qyqLcCrGd1QGGTNSS4E3gn8VFX9\nLUaXBH6Z4dX7CUafr3FTa0yyldE1063dPh9OsnR+9Lka3PIBXM5z7zRdB1w367qWqfmzwBsYvfHX\nd20XMOVO2ozr3Aj8d+Bngc93bYOtGTgX+N9T2gdZM6O1XP8HeCmjJRefB944xHqBi4CvLff/lLE7\nvd32F4BtSx17sDMXll6ANzjdHbFLgPs4uQWGs3QT8G7gxFjbkGt+BfCtJJ9I8sdJPprkrzDQmqvq\nUeA/MAqYx4Cnqmo/A613wmI1bmD0GfyBZT+PQw6XVXOlOclLgE8D76qq74731SjmB/OzJPlF4Imq\negCYept/aDUz+u3/GuDDVfUa4GkmTimGVHOSlzL685iLGH0oX5LkV8fHDKnexZxEjUvWP+RweRTY\nNLa9iecm5yB0t9Q/DfxeVf1grc+xJBd0/S8HnphVfVP8XWBHkm8C/xl4XZLfY9g1HwWOVtVXu+1P\nMQqbxwda8xuAb1bVk1X1DPAHjE7zh1rvuMXeB5Ofx41d26KGHC7PLsBLcjaji0n7ZlzTcyQJcBvw\nUFXdPNb1gwWGsPgCw5moqvdU1aaqegWji4x/WFW/xrBrfhw40i22hNGH9+uMrmUMseY/B7Yl+dHu\nPfIGRhfPh1rvuMXeB/uAX05ydpJXMPrbwK8seaRZX1Ba5mLTFcA3GP1pwPWzrmdKfX+P0XWLPwEe\n6B7zwHmMLpg+DNwDrJt1rYvU/1pgX/d80DUDPwl8FXiQ0Uzg3CHXDOwBDgBfY/StAWcNrV5GM9fH\nGC1iPQJctVSNwHu6z+JB4OeXO76L6CQ1MeTTIkmrmOEiqQnDRVIThoukJgwXSU0YLpKaMFwkNWG4\nSGri/wPAKcxOev7ujwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e2d5c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "viewer = fp.Viewer(vars=C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us solve\n",
    "$$\\frac{\\partial C}{\\partial t} = \\nabla\\cdot\\left[\\left(3 + C^3\\right)\\left(2 + C^2\\right)\\nabla C\\right]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Cface = C.faceValue\n",
    "eq = fp.TransientTerm(var=C) == fp.DiffusionTerm(coeff=(3 + Cface**3) * (2 + Cface**2), var=C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "C.constrain(1., where=mesh.facesLeft)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "C.value = 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fixedPoint = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "C.updateOld()\n",
    "for sweep in xrange(20):\n",
    "    res = eq.sweep(dt=1000.)\n",
    "    fixedPoint.append([sweep, res])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fixedPoint = fp.numerix.array(fixedPoint)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEKCAYAAAAM4tCNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHaFJREFUeJzt3XmUVNW1x/HvphEiyiCioEAEFBEc4oAMapIyorboCw7L\npxhHHIhKgkQTgsmLrTEi+oxxTjQOiQmiEZEhRkFjuRyIwgMBBRQEIoPDgzggioK93x/nQuq1TVd1\ndVffW1W/z1q1rDvWLuzefc65ZzB3R0SkvprFHYCIFCclDxHJi5KHiORFyUNE8qLkISJ5UfIQkbwo\neYhIXpQ8JC9mdoaZzTaz9Wa2xsyeMLPD445Lmo6Sh9Sbmf0IuBm4FtgV6ArcAXw3zrikaZl6mEp9\nmFlbYBVwrrtPjDseiY9KHlJfA4GvAZPiDkTipeQh9bUzsNbdq+MOROKl5CH1tQ7oYGb62Slz+gGQ\n+poJfA6cFHcgEi8lD6kXd/8I+AVwh5kNMbNWZradmR1nZuPijk+ajp62SF7M7AxgFNAbWA/MBn7l\n7v+INTBpMkoeIpIXVVtEJC9KHiKSFyUPEcmLkoeI5KV53AGYmVpsRWLi7pbvtYkoecyZ43Ts6Hz4\noeOe/NdVV10VewylHK9ibppXQyUieRx0EJxwAlx7bdyRiEiuEpE8ICSO+++HpUvjjkREcpGY5NGp\nE1xxRXglXSqVijuEeim2eEExF4PYe5iamW+JYeNG6NMH7rkHjjoq1rBESp6Z4Q1oME1U8gB47DG4\n6ip47bW8v1PixP1vLFKbkkse7qHU8eyzVhK/dNH/oLjDEPmKkkseAPPmwYEHlsYvnZKHJFVJJo9o\nf0n80pXK95DS09DkkZinLSJSXJQ8RCQvSh4NMH78ePr27Uvr1q3ZfffdGTx4MC+++GLcYYk0CSWP\nPP36179m1KhR/PznP+f9999n5cqVXHrppUyZMiXu0ESahBpM8/DRRx/RpUsXHnjgAU455ZQ6z03y\n95DypgbTGMycOZONGzdy0klafUDKV+zzeTSENUIn1HwKBevWraNDhw40a6bcK+WrqJNHXLWBnXfe\nmbVr11JdXa0EImVLP/l5GDhwIC1btmTSJK31LOVLySMPbdu25ZprruHSSy9l8uTJfPrpp2zatIm/\n/e1vjB49Ou7wRJqEnrY0wPjx47n55ptZtGgRrVu3pm/fvvzsZz9jwIABW88phu8h5UljWxKuVL6H\nlB49qhWRWCh5iEhesiYPM6s0s8VmtsTMvtIaaGY7mdkkM5tnZi+b2b4Zx1aY2Xwzm2tmrzR28CIS\nnzrbPMysAngDGASsBmYBQ919UcY5NwIfu/svzawXcIe7D4qOLQcOcfd/1fEZavMQiUGh2zz6AUvd\nfYW7bwImAENqnNMbeBbA3d8AupnZLpkx5huciCRXtuTRGViZsb0q2pdpHnAygJn1A/YAukTHHHja\nzGab2YUND1dEkiJb9/RcytvXA7eY2VxgATAX+DI6doS7r4lKIjPMbLG7P59rcNYYg1dEpCCyJY/V\nQNeM7a6E0sdW7r4eGLZlO2rnWBYdWxP993/NbBKhGvSV5FFVVbX1fSqVIpVK1aud4JNPoHdvmDAB\nDj8858tEyko6nSadTjfa/bI1mDYnNJgeBawBXuGrDaZtgc/c/YuoanK4u59rZq2ACndfb2Y7ANOB\nq919eo3PqLXBtL4eeghuuAFmz4aKigbfTqTkFbTB1N03AyOAp4CFwMPuvsjMhpvZ8Oi0PsACM1sM\nHAuMjPZ3BJ43s1eBl4FpNRNHYzr9dGjTBu6+u1CfICKZEts9PR/z58OgQbBwIXTo0Ci3FClZJTu2\nJV8jR8Jnn6kEIpKNkkcNH34YGk+nToW+fRvttiIlRwPjamjXDq67Di69FKqr445GpHSVXPIAOOcc\naNYM7rsv7khESlfJVVu2mDMHjjsOFi2C9u0b/fYiRU9tHnXYUnW5666C3F6kqCl51OGDD0Lj6V//\nCoccUpCPEClaajCtw047wdixcMklajwVaWwlnTwgNJ5WVMC998YdiUhpKelqyxavvgrHHht6nu68\nc0E/SqRoqM0jRz/8IWzcqJ6nIlsoeeToww+hTx+YNAn69y/4x4kknhpMc9SuXRiyf8kl8OWX2c8X\nkbqVTfIA+N73oHVr+O1v445EpPiVTbVli9dfh1QKFiyATp2a7GNFEkdtHnkYPRrWrIEHH2zSjxVJ\nFCWPPGzYEBpPH3gAjjyyST9aJDHUYJqHHXaAW24JjadffBF3NCLFqSyTB8CQIbDXXnDTTXFHIlKc\nyrLassXy5XDooTBrFnTvHksIIrFRtaUBuneHyy+HESNAy8mK1E9ZJw8IyWP58tDzVERylzV5mFml\nmS02syVmNrqW4zuZ2SQzm2dmL5vZvrlemwQtWoROYyNHwvr1cUcjUjyyrRhXQVgxbhBh6clZfHXF\nuBuBj939l2bWC7jD3Qflcm10fWxtHpnOOy90Yb/55rgjEWkahW7z6AcsdfcV7r4JmAAMqXFOb+BZ\nAHd/A+hmZrvmeG1i3HgjjB8f5j4VkeyyJY/OwMqM7VXRvkzzgJMBzKwfsAfQJcdrE6NDB7j+ehg+\nXAPnRHLRPMvxXOoT1wO3mNlcYAEwF/gyx2sBqKqq2vo+lUqRSqVyvbRRnXtu6HV6553wgx/EEoJI\nwaTTadLpdKPdL1ubxwCgyt0ro+0xQLW7j6vjmuXA/sB+uVyblDaPLRYtgm9+E+bNg86JLSeJNFyh\n2zxmAz3NrJuZtQBOA6bUCKBtdAwzuxB4zt0/yeXaJOrdO3RbHzky7khEkq3O5OHum4ERwFPAQuBh\nd19kZsPNbHh0Wh9ggZktBo4FRtZ1bWG+RuO68spQ8pg6Ne5IRJKrrLun1+WZZ2DYsDD/x447xh2N\nSOPTkPwCOvfcsPaL+n5IKVLyKKC1a2G//WDaNOjbN+5oRBqXBsYVUIcOofPYhRfC5s1xRyOSLEoe\nWZx5Zkgiv/lN3JGIJIuqLTl4662w1ssrr0CPHnFHI9I4VG1pAnvuCT/5CVx8seb9ENlCySNHo0bB\ne+/Bn/4UdyQiyaBqSz3Mng0nnBDWfNlll7ijEWkYPaptYldcAe++qxKIFD8ljya2YQPsvz/cfjsM\nHhx3NCL5U4NpE9thB7j77tB4qmkLpZyp5JGnYcOgVatQAhEpRqq2xOSDD0LX9YcfhiOOiDsakfpT\ntSUmO+0Et94KF1wAGzfGHY1I01PyaIBTToF994Vrrok7EpGmp2pLA737LhxwADz5JBx8cNzRiORO\n1ZaYdeoURt4OGwabNsUdjUjTUfJoBGefDbvtBuO2OS20SOlRtaWRvP02HHIIpNOhHUQk6VRtSYiv\nfx2uvTYsW6mJg6QcKHk0ogsvhNatNeeplAdVWxrZ8uVw6KHwwguwzz5xRyOybQWvtphZpZktNrMl\nZja6luNtzWyqmb1qZq+Z2bkZx1aY2Xwzm2tmr+QbZDHp3j30+zjvPK15K6Ut23KTFcAbwCBgNTAL\nGJq5eJOZXQm0dvcxZtYhOr+ju2+Olp48xN3/VcdnlFTJA6C6Go46Co4/PgzhF0miQpc8+gFL3X2F\nu28CJgBDapxTDbSJ3rcB1kWrxW2NMd/gilWzZnDvvXD99bB4cdzRiBRGtuTRGViZsb0q2pfpdqCP\nma0B5hEtNxlx4Gkzmx2tY1s2evSAq69W9UVKV/Msx3OpT1QCc9z9SDPbE5hhZt9w9/XA4e7+jpnt\nEu1f7O7P17xBVVXV1vepVIpUKpXzF0iyiy+GRx+Fm24KEyiLxCmdTpNOpxvtftnaPAYAVe5eGW2P\nAardfVzGOdOAse7+YrT9DDDa3WfXuNdVwCfuflON/SXX5pFp+XLo1w+eew769Ik7GpF/K3Sbx2yg\np5l1M7MWwGnAlBrnvE1oUMXMOgK9gGVm1srMWkf7dwCOARbkG2ix6t49dB475xyNfZHSUmfyiBo+\nRwBPAQuBh919kZkNN7Ph0Wm/BA4zs/nA08BPoqcrnYDnzexV4GVgmrtPL9QXSbKLLoL27UMDqkip\nUCexJrJyZRiyP306HHRQ3NGIaGxL0ejaNTScnn02fP553NGINJxKHk3IHU4+GXr1UhVG4qcJkIvM\n++/DN74RHuEefnjc0Ug5U7WlyOy6K9x1V3j68skncUcjkj+VPGJy3nnQsiX89rdxRyLlStWWIvXR\nR6H6cuedWrZS4qHkUcSeew7OOAPmzYMOHeKORsqNkkeR+/GP4a23YOJEsLIbfyxxUoNpkbv2Wli6\nFP7wh7gjEakflTwSYP78MHnQyy+HofwiTUEljxJwwAEwZgycdZZmXpfioeSREJddBttvD2PHxh2J\nSG5UbUmQVavCwlFTp4Y5QEQKSdWWEtKlS+j38b3vqfepJJ9KHgl0/vlhEN1998UdiZQylTxK0C23\nhEWjHnkk7khEtk0lj4SaPTt0W581C/bYI+5opBSp5FGi+vYNC0adeaYe30oyKXkk2BVXwNe+Fnqh\niiSNqi0J9847Ye7TRx6Bb34z7miklKjaUuJ22y0sXXnmmfCvba74K9L0VPIoEqNGwYoV8NhjGn0r\njaPgJQ8zqzSzxWa2xMxG13K8rZlNNbNXzew1Mzs312sld9dfD2+/DXfcEXckIkG25SYrgDcIK8Kt\nBmYBQ919UcY5VwKt3X2MmXWIzu9IWOe2zmuj61XyyNHSpTBwIMyYAQceGHc0UuwKXfLoByx19xXu\nvgmYAAypcU410CZ63wZYF600l8u1Ug977RU6kP3nf8L69XFHI+UuW/LoDKzM2F4V7ct0O9DHzNYA\n84CR9bhW6umMM+Bb34Lvfz90YReJS/Msx3P58awE5rj7kWa2JzDDzL5RnyCqqqq2vk+lUqRSqfpc\nXnZuvRX69w9PYS64IO5opFik02nS6XSj3S9bm8cAoMrdK6PtMUC1u4/LOGcaMNbdX4y2nwFGExJT\nnddG+9XmkYdFi0IJ5JlnwmRCIvVV6DaP2UBPM+tmZi2A04ApNc55m9Aoipl1BHoBy3K8VvLUuzf8\n+tdw6qlq/5B4ZO3nYWbHAb8BKoB73X2smQ0HcPffmdluwAPAboARSiHjt3VtLfdXyaMBLrgAPv0U\n/vxn9f+Q+tHSC2Xus89gwAC4+OLQiCqSKyUP4c03w6LZTz4ZpjEUyYXGtgh77x0Wzz71VPjgg7ij\nkXKhkkcJueyysPrc5MnQTH8WJAuVPGSrG24II2+1fIM0BZU8Sszq1XDooWH5yqOPjjsaSTKVPOT/\n6dwZxo+Hs88Oo3BFCkXJowSlUnD55XDKKbBxY9zRSKlStaVEucNpp0GbNvD738cdjSSRqi1SK7Mw\ncG7mTLj77rijkVKkkkeJW7IkdCCbPDlMJCSyhUoeUqeePeH++0MHsnfeiTsaKSVKHmXg+ONh+PDQ\ngPr553FHI6VC1ZYyUV0dSh/t24c2EI3AFVVbJCfNmsEDD8A//hHGwYg0lEoeZeatt0ID6kMPwZFH\nxh2NxEklD6mXPfcMPVCHDoVly+KORoqZkkcZ+s534Oc/h+9+Fz7+OO5opFip2lKm3MPsY6tXw+OP\nQ0VF3BFJU1O1RfJiBrfdBhs2wJgxcUcjxUjJo4xttx385S8waRLcd1/c0Uixybbok5S4nXeGadPC\nGjA9eoQRuSK5UMlD6NUrPIE5/fQwFkYkF1mTh5lVmtliM1tiZqNrOX6Fmc2NXgvMbLOZtYuOrTCz\n+dGxVwrxBaRxHHUUXHtt6Mq+bl3c0UgxyLbcZAXwBmFFuNXALGCouy/axvknAJe5+5YV5JYDh7j7\nv+r4DD1tSZDRo8Mw/hkzoGXLuKORQir005Z+wFJ3X+Hum4AJwJA6zj8DeKhmjPkGJ01v7Fjo2BGG\nDQvjYUS2JVvy6AyszNheFe37CjNrBRwLTMzY7cDTZjbbzC5sSKDSNJo1gz/+EZYvh//6r7ijkSTL\n9rSlPvWJ/wBecPcPM/Yd7u7vmNkuwAwzW+zuz9e8sKqqauv7VCpFSk3+sdp++zB50GGHQbducKHS\nfklIp9Ok0+lGu1+2No8BQJW7V0bbY4Bqdx9Xy7mTgIfdfcI27nUV8Im731Rjv9o8EmrJkvAI9/e/\nDw2pUloK3eYxG+hpZt3MrAVwGjClliDaAt8CJmfsa2VmraP3OwDHAAvyDVSaXs+eoev6uefCrFlx\nRyNJU2fycPfNwAjgKWAhoWSxyMyGm9nwjFNPBJ5y988y9nUEnjezV4GXgWnuPr1xw5dC698/TKQ8\nZAgsXRp3NJIkGhgnObn7bhg3Dl56KTyNkeLX0GqLuqdLTi66KEygPHgwPPtsWA9GyptKHpIzd7jk\nEnjzTfjrX+FrX4s7ImmIhpY8lDykXr78MoyBqa6GRx7RPCDFTPN5SJOqqIA//Qk++igs56C8X76U\nPKTeWrYMc4AsWAA/+YkSSLlS8pC8tG4Nf/tbeF13XdzRSBz0tEXy1r49TJ8eeqG2aQM/+EHcEUlT\nUvKQBtl9d3j66ZBAdtwRzjsv7oikqSh5SIN16xbm/zjyyPD4dujQuCOSpqDkIY2iV69QhTn66NCg\nevLJcUckhabkIY1mv/3giSegshKaNw+LSknp0tMWaVQHHRR6n15wQfivlC4lD2l0ffvC1Kmh8fSJ\nJ+KORgpFyUMKon9/mDIlzAWiEkhpUvKQghkw4N8lkGnT4o5GGpuShxRU//4hcZx/fujSLqVDT1uk\n4Pr1C93YBw+GL76A006LOyJpDEoe0iQOPjj0A6mshE8/VU/UUqDkIU3mgAPCLGRHHw0bNsCIEXFH\nJA2h5CFNqlcveO65kEDWr4ef/hRMawoWJc0kJrF4552QQAYPDhMrK4E0PU1DKEVr3bqQPPbbD373\nu9ClXZpOwachNLNKM1tsZkvMbHQtx68ws7nRa4GZbTazdrlcK+Vt553hmWdg1So45RT47LPs10hy\nZFtusgJ4AxgErAZmAUPdfdE2zj8BuMzdB+V6rUoe8sUXoSfq22+HXqnt28cdUXkodMmjH7DU3Ve4\n+yZgAjCkjvPPAB7K81opUy1ahEmVBwyAI44ISUSSL1vy6AyszNheFe37CjNrBRwLTKzvtSLNmsF/\n/3cYjXv44fDqq3FHJNlka6KqT33iP4AX3P3D+l5bVVW19X0qlSKVStXjY6WU/OhH0LUrHHMMPPgg\nHHts3BGVjnQ6TTqdbrT7ZWvzGABUuXtltD0GqHb3cbWcO4mwEPaE+lyrNg+pzYsvhkbUq68O68NI\n4yvoo1oza05o9DwKWAO8Qu2Nnm2BZUAXd/+sntcqeUitli6F44+HE06AG27Q6nSNraANpu6+GRgB\nPAUsJJQsFpnZcDPL/HtwIvDUlsRR17X5BirlZ6+9YOZMmDMHTjwRPv447ogkkzqJSeJt2hTWhHnh\nhfAot0ePuCMqDVqrVkredtvBXXfB978Phx0WOpZJ/FTykKLy97/DGWeEAXUjR2pMTENobIuUnRUr\n4KSToE8fuOceaNUq7oiKk6otUna6dQuPcisqQq/UJUvijqg8KXlIUWrVCv7wB7j44tAjVfOjNj1V\nW6TozZoFp54aqjLjxoWxMpKdqi1S9g49NPQFWbYsDKxbtizuiMqDkoeUhPbt4fHHw5OYAQNgwoS4\nIyp9qrZIyZkzB4YOhYED4bbboHXruCNKJlVbRGo4+GD4n/8JncsOPBBeeinuiEqTSh5S0iZNCk9k\nhg2Dq66Cli3jjig5VPIQqcNJJ8G8efD662Hlurlz446odCh5SMnr2DE0pl5+eZhc6Be/gM8/jzuq\n4qfkIWXBDM4+O0xvOH8+HHSQ2kIaSm0eUnbc4dFHw8C6E0+E666Ddu3ijqrpqc1DpJ7MQo/U11+H\nL78MA+weeigkFcmdSh5S9l56CS65JHQ0u+022HffuCNqGip5iDTQYYfB7Nlw8smQSoXqzAcfxB1V\n8il5iBDWyR0xAhYuDCvY7bMP3HpreC+1U/IQybDLLmHKw6efhieeCItwP/qo2kNqozYPkTpMnx6m\nPKyoCE9lBg0qnakPNQ2hSIFVV8Mjj4Tu7Z06wTXXwLe/HXdUDVfwBlMzqzSzxWa2xMxGb+OclJnN\nNbPXzCydsX+Fmc2Pjr2Sb5AicWrWDE4/PTzaHTYMzj8/NKw+80x5V2eyrRhXQVj1bRCwGphFjVXf\nzKwd8CJwrLuvMrMO7r42OrYcOMTd/1XHZ6jkIUVl8+bQL+S668Jw/9GjQ2ezYlvRrtAlj37AUndf\n4e6bgAnAkBrnnAFMdPdVAFsSR2aM+QYnkkTNm8NZZ4WSyJVXwo03Qq9ecPvt8MkncUfXdLIlj87A\nyoztVdG+TD2B9mb2rJnNNrOzMo458HS0/8KGhyuSHM2ahRLHzJnwxz/Cs8/CHnvAqFHw5ptxR1d4\nzbMcz6U+sR1wMGFB61bATDP7h7svAY5w9zVmtgsww8wWu/vzNW9QVVW19X0qlSKVSuUYvkj8zEJH\ns8MOg3/+MzzqPeII2H9/uOiikGCSMI9IOp0mnU432v2ytXkMAKrcvTLaHgNUu/u4jHNGA9u7e1W0\n/XvgSXd/tMa9rgI+cfebauxXm4eUnM8/DxMR3X13GMU7dGgY1du3b3Ie9Ra6zWM20NPMuplZC+A0\nYEqNcyYDR5hZhZm1AvoDC82slZm1joLcATgGWJBvoCLFpGXL8ITm738PS0N06BAmZ95nn/DI9/XX\n446w4bL28zCz44DfABXAve4+1syGA7j776JzrgDOA6qBe9z9VjPrATwW3aY58Gd3H1vL/VXykLLg\nDq+8EvqM/OUvsMMOoUozZEhYPqKpn9aok5hIEaquDpM0P/44TJ4M778PlZVwzDGhF2unToWPQaNq\nY9CYjU5NodjihdKPuVmzUNr41a/gtddCiWTgQHjsMejdO8wxcvHFMH48LF+ezM5oSh55KLYf7GKL\nF8ov5m7dQrJ47DFYuxYefBB69oSJE8NavLvtBscfH+ZfnTgxLO795ZeNFnpesj2qFZEmVlEBhxwS\nXj/6USh1rFwZqjlz5oQFvufPh/feg732gr33hh49QgLq2hW6dAmTPu+6a1i7plCUPEQSzgy+/vXw\nOumkf+/fsAHeeAPeeiusz7tgQZhGYPXqkFjWroXtt4eddoI2bWDHHcN2ixaNM1taIhpMYw1ApIwV\n9dMWESlOajAVkbwoeYhIXmJNHrlMNBQnM+sajRZ+PZro6IfR/vZmNsPM3jSz6dGcJokSDReYa2ZT\no+3Exmxm7czsUTNbZGYLzax/kuMFMLNR0c/EAjMbb2Ytkxazmd1nZu+Z2YKMfduM0czGRL+Li83s\nmGz3jy15RBMN3Q5UAn2AoWbWO654tmETMMrd9wUGAJdGMf4UmOHuewPPRNtJMxJYyL9HRic55luA\nJ9y9N3AAsJgEx2tmnYEfECa62p8wdON0khfz/YTfr0y1xmhmfQhj1/pE19xpZnXnB3eP5QUMJIy+\n3bL9U+CnccWTY8yPE2ZVWwx0jPZ1AhbHHVuNOLsATwNHAlOjfYmMGWgLLKtlfyLjjeLpDLwN7ETo\n7jAVODqJMQPdgAXZ/l2BMcDojPOeBAbUde84qy25TDSUGGbWDTgIeJnwj/9edOg9oGNMYW3LzcCP\nCQMVt0hqzN2B/zWz+81sjpndE43CTmq8uPtq4CZCAlkDfOjuM0hwzBm2FePuhN/BLbL+PsaZPIrm\nGbGZ7QhMBEa6+/rMYx7SdGK+i5mdALzv7nPZxhSQCYu5OWEyqTvd/WBgAzWK+wmLFzPbCfgu4a/6\n7sCOZnZm5jlJi7k2OcRYZ/xxJo/VQNeM7a78/8yXCGa2HSFxPOjuj0e73zOzTtHx3YD344qvFocB\n340mn34I+I6ZPUhyY14FrHL3WdH2o4Rk8m5C44VQdV3u7uvcfTNh6omBJDvmLbb1c1Dz97FLtG+b\n4kweuUw0FCszM+BeYKG7/ybj0BTgnOj9OYS2kERw9yvdvau7dyc04v3d3c8ioTG7+7vASjPbO9o1\nCHid0I6QuHgj/wQGmNn20c/IIELjdJJj3mJbPwdTgNPNrIWZdSfMTVz3cikxN+YcR1jaYSkwJu7G\npVriO4LQbvAqMDd6VQLtCQ2SbwLTgXZxx7qN+L8NTIneJzZm4BuEZT3mEf6Kt01yvFHMVcAiwux4\nfyDM5ZuomAklzzXAF4T2xfPqihG4MvpdXExYSqXO+6t7uojkRT1MRSQvSh4ikhclDxHJi5KHiORF\nyUNE8qLkISJ5UfIQkbwoeYhIXv4PO/cWKzNYvu4AAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1107058d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For Newton's method, we solve solve the set of equations $\\mathbf{F}(\\mathbf{x})$ on the variables $\\mathbf{x}$\n",
    "by the first order Taylor expansion\n",
    "$$ \\mathbf{F} + \\mathbf{J}\\cdot\\delta \\mathbf{x} = 0$$\n",
    "where $\\mathbf{J}$ is the Jacobian of $\\mathbf{F}$:\n",
    "$$J_{ij} = \\frac{\\partial F_i}{\\partial x_j}$$\n",
    "and $\\delta \\mathbf{x}$ is the variation in $\\mathbf{x}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The way I get the Jacobian term may not be 100% kosher, but it seems to work. \n",
    "\n",
    "Let us find\n",
    "$$\\mathbf{F}(\\mathbf{x} + \\delta \\mathbf{x}) = \\mathbf{F}(\\mathbf{x}) + \\left.\\delta \\mathbf{F}\\right\\rvert_\\mathbf{x} \\approx 0$$\n",
    "where $\\delta$ is a variation *operator*, such that for\n",
    "$$ \\mathbf{F}(\\mathbf{x}) = \\mathbf{F}(C) \\equiv \\frac{\\partial C}{\\partial t} = \\nabla\\cdot\\left[\\left(3 + C^3\\right)\\left(2 + C^2\\right)\\nabla C\\right]$$\n",
    "then\n",
    "$$\\begin{align*}\n",
    "\\delta \\mathbf{F}(C) &\\equiv \\delta\\left\\{\\frac{\\partial C}{\\partial t} = \\nabla\\cdot\\left[\\left(3 + C^3\\right)\\left(2 + C^2\\right)\\nabla C\\right]\\right\\} \\\\\n",
    "&\\equiv \\frac{\\partial\\, \\delta C}{\\partial t} = \\nabla\\cdot\\left\\{\\delta C\\left[3 C^2 \\left(2 + C^2\\right) + \\left(3 + C^3\\right) 2 C\\right] \\nabla C\\right\\}\n",
    "+ \\nabla\\cdot\\left[\\left(3 + C^3\\right)\\left(2 + C^2\\right)\\nabla \\delta C\\right]\n",
    "\\end{align*}$$\n",
    "\n",
    "I think there's probably a more proper way to get here via variational derivatives and the Euler-Lagrange equation, but it's leaving me with a couple of extra terms that blow up the solution. I'll try to find time to work through it more rigorously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "deltaC = fp.CellVariable(mesh=mesh, name=r\"$\\delta C$\", hasOld=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "newtonEq = ((fp.TransientTerm(var=deltaC) \n",
    "             == fp.ConvectionTerm(coeff=(3*Cface**2*(2 + Cface**2)\n",
    "                                         + (3 + Cface**3)*2*Cface) * C.faceGrad, var=deltaC)\n",
    "             + fp.DiffusionTerm(coeff=(3 + Cface**3) * (2 + Cface**2), var=deltaC)) \n",
    "            + fp.ResidualTerm(equation=eq, underRelaxation=1.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "deltaC.constrain(0., where=mesh.facesLeft)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "C.value = 0\n",
    "deltaC.value = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "newton = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "C.updateOld()\n",
    "deltaC.updateOld\n",
    "for sweep in xrange(20):\n",
    "    res = newtonEq.sweep(dt=1000.)\n",
    "    C.value = C.value + deltaC.value\n",
    "    newton.append([sweep, res, max(abs(deltaC))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "newton = fp.numerix.array(newton)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEKCAYAAAAM4tCNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHaFJREFUeJzt3XmUVNW1x/HvphEiyiCioEAEFBEc4oAMapIyorboCw7L\npxhHHIhKgkQTgsmLrTEi+oxxTjQOiQmiEZEhRkFjuRyIwgMBBRQEIoPDgzggioK93x/nQuq1TVd1\ndVffW1W/z1q1rDvWLuzefc65ZzB3R0SkvprFHYCIFCclDxHJi5KHiORFyUNE8qLkISJ5UfIQkbwo\neYhIXpQ8JC9mdoaZzTaz9Wa2xsyeMLPD445Lmo6Sh9Sbmf0IuBm4FtgV6ArcAXw3zrikaZl6mEp9\nmFlbYBVwrrtPjDseiY9KHlJfA4GvAZPiDkTipeQh9bUzsNbdq+MOROKl5CH1tQ7oYGb62Slz+gGQ\n+poJfA6cFHcgEi8lD6kXd/8I+AVwh5kNMbNWZradmR1nZuPijk+ajp62SF7M7AxgFNAbWA/MBn7l\n7v+INTBpMkoeIpIXVVtEJC9KHiKSFyUPEcmLkoeI5KV53AGYmVpsRWLi7pbvtYkoecyZ43Ts6Hz4\noeOe/NdVV10VewylHK9ibppXQyUieRx0EJxwAlx7bdyRiEiuEpE8ICSO+++HpUvjjkREcpGY5NGp\nE1xxRXglXSqVijuEeim2eEExF4PYe5iamW+JYeNG6NMH7rkHjjoq1rBESp6Z4Q1oME1U8gB47DG4\n6ip47bW8v1PixP1vLFKbkkse7qHU8eyzVhK/dNH/oLjDEPmKkkseAPPmwYEHlsYvnZKHJFVJJo9o\nf0n80pXK95DS09DkkZinLSJSXJQ8RCQvSh4NMH78ePr27Uvr1q3ZfffdGTx4MC+++GLcYYk0CSWP\nPP36179m1KhR/PznP+f9999n5cqVXHrppUyZMiXu0ESahBpM8/DRRx/RpUsXHnjgAU455ZQ6z03y\n95DypgbTGMycOZONGzdy0klafUDKV+zzeTSENUIn1HwKBevWraNDhw40a6bcK+WrqJNHXLWBnXfe\nmbVr11JdXa0EImVLP/l5GDhwIC1btmTSJK31LOVLySMPbdu25ZprruHSSy9l8uTJfPrpp2zatIm/\n/e1vjB49Ou7wRJqEnrY0wPjx47n55ptZtGgRrVu3pm/fvvzsZz9jwIABW88phu8h5UljWxKuVL6H\nlB49qhWRWCh5iEhesiYPM6s0s8VmtsTMvtIaaGY7mdkkM5tnZi+b2b4Zx1aY2Xwzm2tmrzR28CIS\nnzrbPMysAngDGASsBmYBQ919UcY5NwIfu/svzawXcIe7D4qOLQcOcfd/1fEZavMQiUGh2zz6AUvd\nfYW7bwImAENqnNMbeBbA3d8AupnZLpkx5huciCRXtuTRGViZsb0q2pdpHnAygJn1A/YAukTHHHja\nzGab2YUND1dEkiJb9/RcytvXA7eY2VxgATAX+DI6doS7r4lKIjPMbLG7P59rcNYYg1dEpCCyJY/V\nQNeM7a6E0sdW7r4eGLZlO2rnWBYdWxP993/NbBKhGvSV5FFVVbX1fSqVIpVK1aud4JNPoHdvmDAB\nDj8858tEyko6nSadTjfa/bI1mDYnNJgeBawBXuGrDaZtgc/c/YuoanK4u59rZq2ACndfb2Y7ANOB\nq919eo3PqLXBtL4eeghuuAFmz4aKigbfTqTkFbTB1N03AyOAp4CFwMPuvsjMhpvZ8Oi0PsACM1sM\nHAuMjPZ3BJ43s1eBl4FpNRNHYzr9dGjTBu6+u1CfICKZEts9PR/z58OgQbBwIXTo0Ci3FClZJTu2\nJV8jR8Jnn6kEIpKNkkcNH34YGk+nToW+fRvttiIlRwPjamjXDq67Di69FKqr445GpHSVXPIAOOcc\naNYM7rsv7khESlfJVVu2mDMHjjsOFi2C9u0b/fYiRU9tHnXYUnW5666C3F6kqCl51OGDD0Lj6V//\nCoccUpCPEClaajCtw047wdixcMklajwVaWwlnTwgNJ5WVMC998YdiUhpKelqyxavvgrHHht6nu68\nc0E/SqRoqM0jRz/8IWzcqJ6nIlsoeeToww+hTx+YNAn69y/4x4kknhpMc9SuXRiyf8kl8OWX2c8X\nkbqVTfIA+N73oHVr+O1v445EpPiVTbVli9dfh1QKFiyATp2a7GNFEkdtHnkYPRrWrIEHH2zSjxVJ\nFCWPPGzYEBpPH3gAjjyyST9aJDHUYJqHHXaAW24JjadffBF3NCLFqSyTB8CQIbDXXnDTTXFHIlKc\nyrLassXy5XDooTBrFnTvHksIIrFRtaUBuneHyy+HESNAy8mK1E9ZJw8IyWP58tDzVERylzV5mFml\nmS02syVmNrqW4zuZ2SQzm2dmL5vZvrlemwQtWoROYyNHwvr1cUcjUjyyrRhXQVgxbhBh6clZfHXF\nuBuBj939l2bWC7jD3Qflcm10fWxtHpnOOy90Yb/55rgjEWkahW7z6AcsdfcV7r4JmAAMqXFOb+BZ\nAHd/A+hmZrvmeG1i3HgjjB8f5j4VkeyyJY/OwMqM7VXRvkzzgJMBzKwfsAfQJcdrE6NDB7j+ehg+\nXAPnRHLRPMvxXOoT1wO3mNlcYAEwF/gyx2sBqKqq2vo+lUqRSqVyvbRRnXtu6HV6553wgx/EEoJI\nwaTTadLpdKPdL1ubxwCgyt0ro+0xQLW7j6vjmuXA/sB+uVyblDaPLRYtgm9+E+bNg86JLSeJNFyh\n2zxmAz3NrJuZtQBOA6bUCKBtdAwzuxB4zt0/yeXaJOrdO3RbHzky7khEkq3O5OHum4ERwFPAQuBh\nd19kZsPNbHh0Wh9ggZktBo4FRtZ1bWG+RuO68spQ8pg6Ne5IRJKrrLun1+WZZ2DYsDD/x447xh2N\nSOPTkPwCOvfcsPaL+n5IKVLyKKC1a2G//WDaNOjbN+5oRBqXBsYVUIcOofPYhRfC5s1xRyOSLEoe\nWZx5Zkgiv/lN3JGIJIuqLTl4662w1ssrr0CPHnFHI9I4VG1pAnvuCT/5CVx8seb9ENlCySNHo0bB\ne+/Bn/4UdyQiyaBqSz3Mng0nnBDWfNlll7ijEWkYPaptYldcAe++qxKIFD8ljya2YQPsvz/cfjsM\nHhx3NCL5U4NpE9thB7j77tB4qmkLpZyp5JGnYcOgVatQAhEpRqq2xOSDD0LX9YcfhiOOiDsakfpT\ntSUmO+0Et94KF1wAGzfGHY1I01PyaIBTToF994Vrrok7EpGmp2pLA737LhxwADz5JBx8cNzRiORO\n1ZaYdeoURt4OGwabNsUdjUjTUfJoBGefDbvtBuO2OS20SOlRtaWRvP02HHIIpNOhHUQk6VRtSYiv\nfx2uvTYsW6mJg6QcKHk0ogsvhNatNeeplAdVWxrZ8uVw6KHwwguwzz5xRyOybQWvtphZpZktNrMl\nZja6luNtzWyqmb1qZq+Z2bkZx1aY2Xwzm2tmr+QbZDHp3j30+zjvPK15K6Ut23KTFcAbwCBgNTAL\nGJq5eJOZXQm0dvcxZtYhOr+ju2+Olp48xN3/VcdnlFTJA6C6Go46Co4/PgzhF0miQpc8+gFL3X2F\nu28CJgBDapxTDbSJ3rcB1kWrxW2NMd/gilWzZnDvvXD99bB4cdzRiBRGtuTRGViZsb0q2pfpdqCP\nma0B5hEtNxlx4Gkzmx2tY1s2evSAq69W9UVKV/Msx3OpT1QCc9z9SDPbE5hhZt9w9/XA4e7+jpnt\nEu1f7O7P17xBVVXV1vepVIpUKpXzF0iyiy+GRx+Fm24KEyiLxCmdTpNOpxvtftnaPAYAVe5eGW2P\nAardfVzGOdOAse7+YrT9DDDa3WfXuNdVwCfuflON/SXX5pFp+XLo1w+eew769Ik7GpF/K3Sbx2yg\np5l1M7MWwGnAlBrnvE1oUMXMOgK9gGVm1srMWkf7dwCOARbkG2ix6t49dB475xyNfZHSUmfyiBo+\nRwBPAQuBh919kZkNN7Ph0Wm/BA4zs/nA08BPoqcrnYDnzexV4GVgmrtPL9QXSbKLLoL27UMDqkip\nUCexJrJyZRiyP306HHRQ3NGIaGxL0ejaNTScnn02fP553NGINJxKHk3IHU4+GXr1UhVG4qcJkIvM\n++/DN74RHuEefnjc0Ug5U7WlyOy6K9x1V3j68skncUcjkj+VPGJy3nnQsiX89rdxRyLlStWWIvXR\nR6H6cuedWrZS4qHkUcSeew7OOAPmzYMOHeKORsqNkkeR+/GP4a23YOJEsLIbfyxxUoNpkbv2Wli6\nFP7wh7gjEakflTwSYP78MHnQyy+HofwiTUEljxJwwAEwZgycdZZmXpfioeSREJddBttvD2PHxh2J\nSG5UbUmQVavCwlFTp4Y5QEQKSdWWEtKlS+j38b3vqfepJJ9KHgl0/vlhEN1998UdiZQylTxK0C23\nhEWjHnkk7khEtk0lj4SaPTt0W581C/bYI+5opBSp5FGi+vYNC0adeaYe30oyKXkk2BVXwNe+Fnqh\niiSNqi0J9847Ye7TRx6Bb34z7miklKjaUuJ22y0sXXnmmfCvba74K9L0VPIoEqNGwYoV8NhjGn0r\njaPgJQ8zqzSzxWa2xMxG13K8rZlNNbNXzew1Mzs312sld9dfD2+/DXfcEXckIkG25SYrgDcIK8Kt\nBmYBQ919UcY5VwKt3X2MmXWIzu9IWOe2zmuj61XyyNHSpTBwIMyYAQceGHc0UuwKXfLoByx19xXu\nvgmYAAypcU410CZ63wZYF600l8u1Ug977RU6kP3nf8L69XFHI+UuW/LoDKzM2F4V7ct0O9DHzNYA\n84CR9bhW6umMM+Bb34Lvfz90YReJS/Msx3P58awE5rj7kWa2JzDDzL5RnyCqqqq2vk+lUqRSqfpc\nXnZuvRX69w9PYS64IO5opFik02nS6XSj3S9bm8cAoMrdK6PtMUC1u4/LOGcaMNbdX4y2nwFGExJT\nnddG+9XmkYdFi0IJ5JlnwmRCIvVV6DaP2UBPM+tmZi2A04ApNc55m9Aoipl1BHoBy3K8VvLUuzf8\n+tdw6qlq/5B4ZO3nYWbHAb8BKoB73X2smQ0HcPffmdluwAPAboARSiHjt3VtLfdXyaMBLrgAPv0U\n/vxn9f+Q+tHSC2Xus89gwAC4+OLQiCqSKyUP4c03w6LZTz4ZpjEUyYXGtgh77x0Wzz71VPjgg7ij\nkXKhkkcJueyysPrc5MnQTH8WJAuVPGSrG24II2+1fIM0BZU8Sszq1XDooWH5yqOPjjsaSTKVPOT/\n6dwZxo+Hs88Oo3BFCkXJowSlUnD55XDKKbBxY9zRSKlStaVEucNpp0GbNvD738cdjSSRqi1SK7Mw\ncG7mTLj77rijkVKkkkeJW7IkdCCbPDlMJCSyhUoeUqeePeH++0MHsnfeiTsaKSVKHmXg+ONh+PDQ\ngPr553FHI6VC1ZYyUV0dSh/t24c2EI3AFVVbJCfNmsEDD8A//hHGwYg0lEoeZeatt0ID6kMPwZFH\nxh2NxEklD6mXPfcMPVCHDoVly+KORoqZkkcZ+s534Oc/h+9+Fz7+OO5opFip2lKm3MPsY6tXw+OP\nQ0VF3BFJU1O1RfJiBrfdBhs2wJgxcUcjxUjJo4xttx385S8waRLcd1/c0Uixybbok5S4nXeGadPC\nGjA9eoQRuSK5UMlD6NUrPIE5/fQwFkYkF1mTh5lVmtliM1tiZqNrOX6Fmc2NXgvMbLOZtYuOrTCz\n+dGxVwrxBaRxHHUUXHtt6Mq+bl3c0UgxyLbcZAXwBmFFuNXALGCouy/axvknAJe5+5YV5JYDh7j7\nv+r4DD1tSZDRo8Mw/hkzoGXLuKORQir005Z+wFJ3X+Hum4AJwJA6zj8DeKhmjPkGJ01v7Fjo2BGG\nDQvjYUS2JVvy6AyszNheFe37CjNrBRwLTMzY7cDTZjbbzC5sSKDSNJo1gz/+EZYvh//6r7ijkSTL\n9rSlPvWJ/wBecPcPM/Yd7u7vmNkuwAwzW+zuz9e8sKqqauv7VCpFSk3+sdp++zB50GGHQbducKHS\nfklIp9Ok0+lGu1+2No8BQJW7V0bbY4Bqdx9Xy7mTgIfdfcI27nUV8Im731Rjv9o8EmrJkvAI9/e/\nDw2pUloK3eYxG+hpZt3MrAVwGjClliDaAt8CJmfsa2VmraP3OwDHAAvyDVSaXs+eoev6uefCrFlx\nRyNJU2fycPfNwAjgKWAhoWSxyMyGm9nwjFNPBJ5y988y9nUEnjezV4GXgWnuPr1xw5dC698/TKQ8\nZAgsXRp3NJIkGhgnObn7bhg3Dl56KTyNkeLX0GqLuqdLTi66KEygPHgwPPtsWA9GyptKHpIzd7jk\nEnjzTfjrX+FrX4s7ImmIhpY8lDykXr78MoyBqa6GRx7RPCDFTPN5SJOqqIA//Qk++igs56C8X76U\nPKTeWrYMc4AsWAA/+YkSSLlS8pC8tG4Nf/tbeF13XdzRSBz0tEXy1r49TJ8eeqG2aQM/+EHcEUlT\nUvKQBtl9d3j66ZBAdtwRzjsv7oikqSh5SIN16xbm/zjyyPD4dujQuCOSpqDkIY2iV69QhTn66NCg\nevLJcUckhabkIY1mv/3giSegshKaNw+LSknp0tMWaVQHHRR6n15wQfivlC4lD2l0ffvC1Kmh8fSJ\nJ+KORgpFyUMKon9/mDIlzAWiEkhpUvKQghkw4N8lkGnT4o5GGpuShxRU//4hcZx/fujSLqVDT1uk\n4Pr1C93YBw+GL76A006LOyJpDEoe0iQOPjj0A6mshE8/VU/UUqDkIU3mgAPCLGRHHw0bNsCIEXFH\nJA2h5CFNqlcveO65kEDWr4ef/hRMawoWJc0kJrF4552QQAYPDhMrK4E0PU1DKEVr3bqQPPbbD373\nu9ClXZpOwachNLNKM1tsZkvMbHQtx68ws7nRa4GZbTazdrlcK+Vt553hmWdg1So45RT47LPs10hy\nZFtusgJ4AxgErAZmAUPdfdE2zj8BuMzdB+V6rUoe8sUXoSfq22+HXqnt28cdUXkodMmjH7DU3Ve4\n+yZgAjCkjvPPAB7K81opUy1ahEmVBwyAI44ISUSSL1vy6AyszNheFe37CjNrBRwLTKzvtSLNmsF/\n/3cYjXv44fDqq3FHJNlka6KqT33iP4AX3P3D+l5bVVW19X0qlSKVStXjY6WU/OhH0LUrHHMMPPgg\nHHts3BGVjnQ6TTqdbrT7ZWvzGABUuXtltD0GqHb3cbWcO4mwEPaE+lyrNg+pzYsvhkbUq68O68NI\n4yvoo1oza05o9DwKWAO8Qu2Nnm2BZUAXd/+sntcqeUitli6F44+HE06AG27Q6nSNraANpu6+GRgB\nPAUsJJQsFpnZcDPL/HtwIvDUlsRR17X5BirlZ6+9YOZMmDMHTjwRPv447ogkkzqJSeJt2hTWhHnh\nhfAot0ePuCMqDVqrVkredtvBXXfB978Phx0WOpZJ/FTykKLy97/DGWeEAXUjR2pMTENobIuUnRUr\n4KSToE8fuOceaNUq7oiKk6otUna6dQuPcisqQq/UJUvijqg8KXlIUWrVCv7wB7j44tAjVfOjNj1V\nW6TozZoFp54aqjLjxoWxMpKdqi1S9g49NPQFWbYsDKxbtizuiMqDkoeUhPbt4fHHw5OYAQNgwoS4\nIyp9qrZIyZkzB4YOhYED4bbboHXruCNKJlVbRGo4+GD4n/8JncsOPBBeeinuiEqTSh5S0iZNCk9k\nhg2Dq66Cli3jjig5VPIQqcNJJ8G8efD662Hlurlz446odCh5SMnr2DE0pl5+eZhc6Be/gM8/jzuq\n4qfkIWXBDM4+O0xvOH8+HHSQ2kIaSm0eUnbc4dFHw8C6E0+E666Ddu3ijqrpqc1DpJ7MQo/U11+H\nL78MA+weeigkFcmdSh5S9l56CS65JHQ0u+022HffuCNqGip5iDTQYYfB7Nlw8smQSoXqzAcfxB1V\n8il5iBDWyR0xAhYuDCvY7bMP3HpreC+1U/IQybDLLmHKw6efhieeCItwP/qo2kNqozYPkTpMnx6m\nPKyoCE9lBg0qnakPNQ2hSIFVV8Mjj4Tu7Z06wTXXwLe/HXdUDVfwBlMzqzSzxWa2xMxGb+OclJnN\nNbPXzCydsX+Fmc2Pjr2Sb5AicWrWDE4/PTzaHTYMzj8/NKw+80x5V2eyrRhXQVj1bRCwGphFjVXf\nzKwd8CJwrLuvMrMO7r42OrYcOMTd/1XHZ6jkIUVl8+bQL+S668Jw/9GjQ2ezYlvRrtAlj37AUndf\n4e6bgAnAkBrnnAFMdPdVAFsSR2aM+QYnkkTNm8NZZ4WSyJVXwo03Qq9ecPvt8MkncUfXdLIlj87A\nyoztVdG+TD2B9mb2rJnNNrOzMo458HS0/8KGhyuSHM2ahRLHzJnwxz/Cs8/CHnvAqFHw5ptxR1d4\nzbMcz6U+sR1wMGFB61bATDP7h7svAY5w9zVmtgsww8wWu/vzNW9QVVW19X0qlSKVSuUYvkj8zEJH\ns8MOg3/+MzzqPeII2H9/uOiikGCSMI9IOp0mnU432v2ytXkMAKrcvTLaHgNUu/u4jHNGA9u7e1W0\n/XvgSXd/tMa9rgI+cfebauxXm4eUnM8/DxMR3X13GMU7dGgY1du3b3Ie9Ra6zWM20NPMuplZC+A0\nYEqNcyYDR5hZhZm1AvoDC82slZm1joLcATgGWJBvoCLFpGXL8ITm738PS0N06BAmZ95nn/DI9/XX\n446w4bL28zCz44DfABXAve4+1syGA7j776JzrgDOA6qBe9z9VjPrATwW3aY58Gd3H1vL/VXykLLg\nDq+8EvqM/OUvsMMOoUozZEhYPqKpn9aok5hIEaquDpM0P/44TJ4M778PlZVwzDGhF2unToWPQaNq\nY9CYjU5NodjihdKPuVmzUNr41a/gtddCiWTgQHjsMejdO8wxcvHFMH48LF+ezM5oSh55KLYf7GKL\nF8ov5m7dQrJ47DFYuxYefBB69oSJE8NavLvtBscfH+ZfnTgxLO795ZeNFnpesj2qFZEmVlEBhxwS\nXj/6USh1rFwZqjlz5oQFvufPh/feg732gr33hh49QgLq2hW6dAmTPu+6a1i7plCUPEQSzgy+/vXw\nOumkf+/fsAHeeAPeeiusz7tgQZhGYPXqkFjWroXtt4eddoI2bWDHHcN2ixaNM1taIhpMYw1ApIwV\n9dMWESlOajAVkbwoeYhIXmJNHrlMNBQnM+sajRZ+PZro6IfR/vZmNsPM3jSz6dGcJokSDReYa2ZT\no+3Exmxm7czsUTNbZGYLzax/kuMFMLNR0c/EAjMbb2Ytkxazmd1nZu+Z2YKMfduM0czGRL+Li83s\nmGz3jy15RBMN3Q5UAn2AoWbWO654tmETMMrd9wUGAJdGMf4UmOHuewPPRNtJMxJYyL9HRic55luA\nJ9y9N3AAsJgEx2tmnYEfECa62p8wdON0khfz/YTfr0y1xmhmfQhj1/pE19xpZnXnB3eP5QUMJIy+\n3bL9U+CnccWTY8yPE2ZVWwx0jPZ1AhbHHVuNOLsATwNHAlOjfYmMGWgLLKtlfyLjjeLpDLwN7ETo\n7jAVODqJMQPdgAXZ/l2BMcDojPOeBAbUde84qy25TDSUGGbWDTgIeJnwj/9edOg9oGNMYW3LzcCP\nCQMVt0hqzN2B/zWz+81sjpndE43CTmq8uPtq4CZCAlkDfOjuM0hwzBm2FePuhN/BLbL+PsaZPIrm\nGbGZ7QhMBEa6+/rMYx7SdGK+i5mdALzv7nPZxhSQCYu5OWEyqTvd/WBgAzWK+wmLFzPbCfgu4a/6\n7sCOZnZm5jlJi7k2OcRYZ/xxJo/VQNeM7a78/8yXCGa2HSFxPOjuj0e73zOzTtHx3YD344qvFocB\n340mn34I+I6ZPUhyY14FrHL3WdH2o4Rk8m5C44VQdV3u7uvcfTNh6omBJDvmLbb1c1Dz97FLtG+b\n4kweuUw0FCszM+BeYKG7/ybj0BTgnOj9OYS2kERw9yvdvau7dyc04v3d3c8ioTG7+7vASjPbO9o1\nCHid0I6QuHgj/wQGmNn20c/IIELjdJJj3mJbPwdTgNPNrIWZdSfMTVz3cikxN+YcR1jaYSkwJu7G\npVriO4LQbvAqMDd6VQLtCQ2SbwLTgXZxx7qN+L8NTIneJzZm4BuEZT3mEf6Kt01yvFHMVcAiwux4\nfyDM5ZuomAklzzXAF4T2xfPqihG4MvpdXExYSqXO+6t7uojkRT1MRSQvSh4ikhclDxHJi5KHiORF\nyUNE8qLkISJ5UfIQkbwoeYhIXv4PO/cWKzNYvu4AAAAASUVORK5CYII=\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110cd7dd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x111031750>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucjHX/x/HXZ3et5JC1lIhGoptORA6xGrfKSlJRUkqK\njkjdRaG6K4eOuqMzcsqh0v0rlNStlo0kOaRbjjVyvp1KCHv4/P64xra2XXZ2Z+aamf08H495mOua\nmet6z5jdz36v73V9v6KqGGOMMYGIczuAMcaY6GPFwxhjTMCseBhjjAmYFQ9jjDEBs+JhjDEmYFY8\njDHGBCzB7QDBJiJ27rExxhSBqkphnxuTLQ9VtVuQbk888YTrGWLlZp+lfZ6RfAtUTBYPY4wxoWXF\nwxhjTMCseJjj8nq9bkeIGfZZBpd9nu6SohzrimQiorH2nowxJtREBA2gwzzmzrYyxoSfSKF/55gI\nEIw/sK14GGOCwlr80SFYhT6qioeIdATaAxWAsar6ucuRjDGmRIrKPg8RqQi8oKo983nM+jyMCTP/\n8XK3Y5hCKOj/KtA+D9fPthKRt0Vkh4iszLM+VURWi8g6ERmQ52WDgVeCnWXVzlXc9MFNrNi+Itib\nNsaYmOJ68QDGAam5V4hIPE5xSAXqA11FpJ44ngVmq+ryYAXYf2Q//T/vz6XjLyXppCTaTW7Hut3r\ngrV5Y4yL1qxZQ4MGDahQoQKjRo3innvuYciQIUHfj8fjYe7cuUHf7pVXXsmkSZOCvt3icr3PQ1XT\nRcSTZ3UTYL2q+gBEZBrQEbgMaANUEJGzVfXNYu6bD378gAfnPIjX4+WHe37gtHKn0aBqA6545wq+\n6vEV1StUL84ujDEue+6552jTpg3Llwft7818iUhIzjr75JNPCv1cr9fLLbfcwh133BH0HHm5XjwK\nUB3YlGt5M9BUVfsAo070Yq/Xi8fjwePx4PV6872YaO3utfSZ3Yct+7bwznXv0OrMVjmP9WrUi72H\n9nL5pMuZ32M+lU+uXOw3ZIxxx8aNG7nkkkvcjhEWgRSvtLQ00tLS8Pl8+Hy+wHfm9mBc/o4bD7Ay\n13InYHSu5W7AqEJuS4/nwJEDOmjuIE1+NllfXPiiHsk8UuBzB3w+QC9+62Ldd2jfcbdpTEl3op87\nt7Ru3Vrj4+P1pJNO0vLly+vatWu1e/fuOnjwYFVVfeaZZ7Rp06aamZmpqqqvvfaannvuuXr48GHN\nysrS4cOHa+3atTU5OVlvuOEG3bNnT862J06cqDVr1tTk5GQdOnSoejwenTt3br45unfvrnfddZde\nfvnlWr58eb300kt148aNOY8vWLBAGzdurKeccopefPHFunDhwpzHLr30Uh0zZoyqqo4bN05btGih\nDz30kCYlJWmtWrV09uzZqqo6cODAnPdarlw57dOnT75ZCvq/8q8v9O/tSOjzyM8WoEau5Ro4rY9i\nmbFmBue+di7r96xnxd0reLD5g5SKL1Xg84e3GU7Dqg3pOK0jhzIPFXf3xpgw++KLL0hJSeHVV19l\n37591KlT55jDS/3796d06dIMGTKEdevWMWjQICZPnkxiYiKjRo1ixowZzJ8/n23btpGUlMR9990H\nwKpVq7j33nuZPHkyW7duZffu3WzefPxfUVOmTOHxxx9n165dNGjQgJtvvhmAPXv20L59e/r168ee\nPXt48MEHad++PXv37gX+ejhs8eLF/O1vf2P37t30798/5xDV0KFDc97r77//zsiRI4P+eR4jkEoT\nqht/bXkkABv86xOB5UC9Qm7rLxV1w54NetWUq/ScUefo5xs+z7fqFiQzK1NveP8G7Ti1o2ZkZQT0\nWmNKivx+7o59PDi3ovB6vTl/uauq3nbbbTktD1VVn8+nlSpV0nr16ukzzzyTs75evXrHtCS2bt2q\npUqV0szMTH3yySe1a9euOY8dOHBAExMTj9vyyP38/fv3a3x8vG7atEknTpyoTZs2Peb5zZs31/Hj\nx+fkHzt2rKo6LY+zzz77mP2KiO7YsSPf95qfgv6viLaWh4hMBRYCdUVkk4j0UNVMoDcwB1gFvKuq\nPwa67UOZh3hq3lM0Gd2EFjVa8P0933PZWZcFtI34uHgmXTuJw1mH6TmjJ9maHWgMY0q8YJWPojpe\nX8CZZ56J1+tl48aNOS0LAJ/Px7XXXktSUhJJSUnUr1+fhIQEduzYwbZt2zjjjDNynnvyySeTnJx8\n3P3nfn7ZsmWpVKkSW7duZdu2bdSsWfMvmbZu3ZrvtqpWrXrMfgH2799fqPcaTK4XD1XtqqrVVLW0\nqtZQ1XH+9bNV9RxVPVtVhwe63U/Xf8r5r5/P8u3LWXrXUh5p+QiJ8YlFypgYn8j066ezbs86Hpzz\n4NEWjjEmBnz88ccsWrSINm3a8NBDD+Wsr1mzJp9++il79+7NuR08eJBq1apx+umns2nTn+f0HDx4\nkN27dxe4D1U95vn79+9nz549VK9enWrVqrFx48Zjnr9x40aqVw/8TM9wjjHmevEIhU7vdaL3J70Z\nmTqSf3f5NzVPqXniF51A2cSyzOo6iy99XzJkfvDPETfGhE7uP/hy39+1axe9evVi7NixjB8/npkz\nZzJ79mwA7r77bgYOHMgvv/wCwM6dO5kxYwYAnTt3ZtasWSxYsIAjR47w+OOPk519/KMSn3zySc7z\nH3vsMZo3b0716tVp164da9euZerUqWRmZvLuu++yevVqrrrqqoDf52mnncaGDRsCfl1RxGTxuPC0\nC/nh3h9oV6ddULebVCaJz7p9xsTvJ/LK4qBf4G6MCZHcf5Hn7oC+6667uOaaa0hNTaVSpUqMHTuW\nnj17snfvXu6//36uvvpqrrjiCipUqEDz5s1ZvHgxAPXr1+fVV1/lpptuolq1alSqVIkaNWrku++j\n+7zpppt48sknSU5OZtmyZbzzzjsAJCcnM2vWLF588UUqV67MCy+8wKxZs6hUqVK+28nbusi9fP/9\n9zN9+nQqVapEv379iv6BFUJUjm11PCKi2dlKKFtvvl99pIxLYXib4XS7oFvodmRMlLCxrY6vR48e\nnHHGGTz99NNuR4mdsa1Cwd/qDBlPRQ9zus3hoc8eYuaamaHdmTEm6sViYY3J4vHkk8U7M6Mw6lep\nz8yuM7ljxh2k+dJCuzNjTFQL1dAlborJw1b16ysvvgipqSd+fnF9+fOXdJnehdk3z6ZRtUah36Ex\nEcgOW0UPO2x1HI8/Hp7WB0DrWq0Z3WE0V029itW7Vod+h8YYEwFisnh07gy//Qafh2mewY5/68gz\nbZ6h7Ttt+eW3X8KzU2OMcVFMFo/4eHjssfC1PgC6N+jOg80e5PJJl/O/A/8r8nZUlYysDLuS3RgT\n0WKyz0NVycqC886DkSPh8svDt/9/pv2TMUvHUL1CdTKzMwO+ZWs2CXEJlI4vTYOqDWh0eiMaV2tM\no2qNOCf5HOLj4sP3ZowpJOvziB7B6vOI2eIBMHkyvP46pKcT0us+clNVlm1fRkZWBglxCQHf4iQO\nEeHXQ7+ybNsylmxdwpJtS/hu63fsOLCDhlUbHlNQ6ibXJU5isgFpoogVj+hhxaMAuYtHVhbUrw+v\nvQZt2rgcLAj2/rGXpduWsmTrEr7b9h1Lti5h18FdXHT6RTkFpXG1xtSuVNsKigkrKx7Rw4pHAXIX\nD4BJk2D0aJg3L3ytj3DafXD3XwrKr4d+5aLTL6LnRT3pel7XmDu/3ESeklo84uLiWL9+PWeddZbb\nUQrNikcB8haPzEyoVw/eegtat3YxWBjtOriLrzd9zeAvB3N6udN546o38FT0uB3LxLCSXDzWrVtH\n7dq13Y5SaHadRyElJMDgwc6ZVyVF5ZMr0+GcDizptQSvx0vjtxoz4usRZGZnuh3NmLDyeDy8+OKL\nXHjhhVSsWJEbb7yRw4cPAzBr1iwaNGhAUlISLVq0YOXKlQCMGzeOq6++OmcbderU4YYbbshZrlGj\nBitWrKBVq1YAXHjhhZQvX573338fgNGjR1OnTh2Sk5Pp2LEj27Zty3ltXFwcb775JnXr1iUpKYne\nvXuH/DMImUBmjnLzBpQFJgBvATcd53maV0aGau3aqmlpf3moRFi7a622Ht9aG73ZSJdtW+Z2HBOD\n8vu5iwQej0ebNm2q27Zt0z179mi9evX0jTfe0KVLl+qpp56qixcv1uzsbJ0wYYJ6PB49cuSIbtiw\nQStWrKiqqlu2bNEzzzxTa9SooaqqGzZs0KSkpJzti4hu2LAhZ3nu3LlauXJlXbZsmR4+fFj79Omj\nrVq1Oub5HTp00N9++01/+eUXrVKlin766adh+jQcBf1fEeBMggluFa0iuA54T1U/FpFpwJTCvjB3\n6+OLL0IXMFLVSa7D3FvnMn75eK6YdAU9GvTgCe8TnFzqZLejmRJCngxOv5s+Efihsb59++bMvteh\nQweWL1/OihUruOuuu7j44osBuPXWWxk2bBiLFi0iJSWF8uXLs2zZMtasWUPbtm1ZsWIFa9asYeHC\nhTktjvxMnjyZO+64gwYNGgAwfPhwkpKS+OWXX3JmC3zkkUeoUKECFSpUoHXr1ixfvpy2bdsG/L7c\n5mrxEJG3gfbA/1T1/FzrU4F/AfHAGFV9FqgOrPA/JSvQfXXrBk8/7Zy2m5JS/OzRRkTo0bAHV9a5\nkn5z+nHB6xfwxlVvBDwtrzFFUZRf+sGSd9rWrVu3smfPHiZMmMCoUaNyHsvIyMiZ+vXSSy8lLS2N\n9evXc+mll1KxYkXmzZvH119/zaWXXlrgvrZt20bjxo1zlsuWLUtycjJbtmzJKR558+SeQjaauN3n\nMQ44ZvhCEYkHXvGvrw90FZF6wGbg6GwrAedOSIBBg0pW30d+Tit3GlM7TeXl1Je5Y8Yd3Pbhbew+\nWPD0mcbEoho1ajBo0KBjppjdv38/Xbp0AZzi8eWXX5Keno7X680pJvPmzTtu8ahWrRo+ny9n+cCB\nA+zevbtIU8pGOleLh6qmA3vzrG4CrFdVn6pmANOAjsC/gU4i8howoyj7u+UW+Okn+Oqr4qSODe3r\ntueHe37glNKncN7r5zFl5ZQSebaMKVmOfsd79erFG2+8weLFi1FVDhw4wMcff5zTCjhaPA4dOkS1\natVo2bIln376KXv27KFhw4Y528s77WvXrl0ZN24cK1as4PDhwwwcOJBmzZrltDoKyhONIrHPozqw\nKdfyZqCpqh4Ebi/MBrxeLx6PB4/Hg9frxev1AlCqFAwc6LQ+wjVoYiQrX7o8L7d7mZvOv4leM3sx\n6ftJvN7+dTut18Sso/NqNGrUiNGjR9O7d2/WrVtHmTJlSElJyWlV1KlTh/Lly5PiP8ZdoUIFateu\nzamnnnrMdVP//Oc/6d69O3/88QejR4+mc+fOPP3003Tq1Im9e/fSokULpk2bdsz+88vjhrS0NNLS\n0vD5fMe0lgrL9es8RMQDzDza5yEinYBUVe3lX+6GUzz6FHJ7erz3dOQI1K0LU6bAJZcUN33syMjK\n4PmFzzPi6xEMShlE36Z9bRwtU2gl9TqPaBTL13ls4c++Dfz3Nwdr44mJf7Y+zJ9KxZdiYMpAFt6x\nkBlrZ9BsbDNWbF9x4hcaY0qkSCweS4A6IuIRkUSgC0Xs4yjIbbfB6tWwaFEwtxob6ibX5Ytbv+Ce\nxvdw2aTLSN+Y7nYkY0wEcrV4iMhUYCFQV0Q2iUgPVc0EegNzgFXAu6r6YzD3a62P4xMRbm94O2Ov\nHkv3D7vz++Hf3Y5kjIkwrvd5BNuJ+jyOOnIEzj4bpk+HJk3CECxK9ZzRE4AxV49xOYmJZNbnET1i\nuc8jLBIT4dFHrfVxIi+1fYkvfv6Cj1Z/5HYUY0wEKbHFA+D22+H77+Hbb91OErnKly7PxGsncvfH\ndxdrel1jTGwpsYetjnr1VZg9G2bNCmGoGPDIfx5h9a7V/F+X/7P5Qcxf2HciugTjsFWJLx6HDjl9\nHx9+CLmGpDF5HM48TJMxTejXtB89GvZwO44xJsiseARYPABGjXKuOJ8R1BOCY8/KHSv5+8S/s7jn\nYmol1XI7jjEmiKx4FKF4HDoEtWvDzJlw0UUhChYjXlj4AjPWzODL7l/aFejGxBA726oITjoJ+veH\np55yO0nke6DZA4gII74e4XYUY4yLrOXh98cfTuvjk0/AP4+LKYDvVx8Xj76YubfO5YLTLnA7jjEm\nCKzlUURlyljro7A8FT08f/nzdPt3Nw5nHnY7jjHGBdbyyOVo62PaNDjOTJMG51S/6967jrqV6vLs\n5c+6HccYU0zW8iiGMmVg7Fi4+WbYudPtNJFNRHjrqreY9P0kGzzRmBLIikce7do5853fcgtkZ7ud\nJrJVKVuFtzq8xa0f3sq+w/vcjmOMCSM7bJWPzEzweqF9e2f8K3N8vWb0IluzGdtxrNtRjDFFZIet\ngiAhwen3ePllSLcjMic0ou0I0jam2eCJxpQgVjwKcMYZMG4c3HST9X+cSPnS5Zl4jQ2eaExJElWH\nrUSkI9AeqACMVdXP83lOsQ9b5fboo7BsmXP9R5yV2uN69D+PsmrXKj7s8qENlGdMlCkRw5OISEXg\nBVXtmc9jQS0e1v9ReEeyjtBkdBP6Nu3L7Q1vdzuOMSYAUdHnISJvi8gOEVmZZ32qiKwWkXUiMuA4\nmxgMvBLalA7r/yi8xPhE3rnuHQb8ZwA/7/3Z7TjGmBBy60DMOCA19woRiccpCKlAfaCriNQTkVtE\n5CURqSaOZ4HZqro8XGFz93/8zw7pH9d5p57HIy0eofuH3cnKznI7jjEmRFwpHqqaDuzNs7oJsF5V\nfaqaAUwDOqrqJFV9QFW3An2ANkBnEbkrnJnt+o/Ce6D5A8THxfPi1y+6HcUYEyIJbgfIpTqwKdfy\nZqBp7ieo6khg5Ik25PV68Xg8eDwevF4vXq83KAGfftrp/3jmGRg4MCibjElxEsf4juO5ePTFpJ6d\naoMnGhOB0tLSSEtLw+fz4fP5An69ax3mIuIBZqrq+f7lTkCqqvbyL3cDmqpqnwC3G9QO87w2b3Zm\nHHzvPRv/6kQmLJ/Ay9+8zHd3fmdnXxkT4aKiw7wAW4AauZZr4LQ+Ior1fxTerRfeiqLMXj/b7SjG\nmCCLpOKxBKgjIh4RSQS6ABE5MWy7dk7fh/V/HJ+I0P+S/jy34Dm3oxhjgsytU3WnAguBuiKySUR6\nqGom0BuYA6wC3lXVH93IVxhPPw0HDzr9H6Zg1597Pb5ffXyz+Ru3oxhjgigqLxI8nlD3eeRm/R+F\nM+qbUczbOI/pN0x3O4oxpgDR3OcRdc44A8aPt/6PE7m94e3M3zifdbvXuR3FGBMkVjyKKTUVbr3V\n+j+Op2xiWe5pfI9d92FMDLHDVkGQmQmtWzsd6Xb9R/52HtjJOa+cw4/3/chp5U5zO44xJg87bOWC\nhASYOhVGjoT5891OE5mqlK1C1/O6MmrxKLejGGOCwFoeQfTpp9CzJyxdCqee6kqEiLZhzwaajW3G\nz/f/TLnEcm7HMcbkYi0PFx3t/+jaFQ4ccDtN5KldqTZ/r/V3xiwd43YUY0wxWfEIsqeegurVISXF\nOZXXHOvhSx5mxNcjyMjKcDuKMaYYrHgEWUICTJgAXbpA06bwjV0bd4zG1RpTJ7kO7/73XbejGGOK\nwYpHCIjAgAHw+utw1VUwZYrbiSLLw5c8zHMLniPW+tuMKUmseITQ1VfD3LnO6buDB9t1IEe1rd0W\ngDkb5ricxBhTVFY8QuyCC2DxYkhLg86drSMd/AMmtrABE42JZlY8wuDUU50WSIUK0LIlbNp04tfE\nui7ndmH9nvUs2brE7SjGmCKw4hEmpUs784DcfDM0awaLFrmdyF2l4kvxYPMHeX7h825HMcYUgV0k\n6IJZs+D222HECGde9JJq/5H91Hq5FovuWETtSrXdjmNMiRboRYJWPFzyww9Oh/qNN8KQIRBXQtuA\ng78YzN4/9vJq+1fdjmJMiWbFI0qKB8DOndCpEyQnw6RJUK4EjtixY/8O6r1ajzW911ClbBW34xhT\nYsX88CQiUlZEvhWR9m5nKa4qVeA//4FKlaBFC9i40e1E4XdaudO44dwbeGXxK25HMcYEIOqKB9Af\niJnLkxMTYcwY6N4dmjeHr792O1H4/aP5P3h9yescOGLnMRsTLdyaw/xtEdkhIivzrE8VkdUisk5E\nBuTzustx5jffGa6s4SACDz7oFJGOHZ1DWCVJneQ6tDqzFW8ve9vtKMaYQnKlz0NEUoD9wERVPd+/\nLh5YA1wGbAG+BboCjYGLgOeBe4GyQH3gD+DavB0c0dTnkZ9Vq6BDB+eCwmHDID7e7UTh8c3mb+gy\nvQvr+64nIS7B7TjGlDhR0eehqunA3jyrmwDrVdWnqhnANKCjqk5S1QdUdauqDlbVB4ApwFtRXSUK\nUL++M5ji4sVOK+S339xOFB5Nz2iKp6KH9//7vttRjDGFEEl/4lUHcl97vRlomt8TVXXC8Tbk9Xrx\neDx4PB68Xi9erzd4KcOgcmX47DPo18+5oHDGDKhTx+1Uode/RX8GfTGIG8+7EZFC/wFkjCmCtLQ0\n0tLS8Pl8+Hy+gF/v2qm6IuIBZuY6bNUJSFXVXv7lbkBTVe0T4HZjqkHy5pvw+OPwzjtw+eVupwkt\nVeWCNy5gxBUjuLx2jL9ZYyJM0A5bich+Efm9gNu+4MQ9xhagRq7lGjitjxLtrrvg/fedGQpffhli\nqC7+hYg4w7UvtAETjYl0BRYPVS2nquULuFUIQZYlQB0R8YhIItAFmBGC/USdVq2cU3jffhvuuAMO\nH3Y7UejceN6NrN61mqXblrodxRhzHIXuMBeRU0Wk5tFbcXYqIlOBhUBdEdkkIj1UNRPoDczBOR33\nXVX9sTj7iSUeDyxYAPv2QevWsH2724lCIzE+kQeaPWADJhoT4U7Y5yEiVwMvAtWA/wFnAj+q6rmh\njxe4WOvzyCs7G55+GsaOhf/7P2jUyO1Ewff74d+p9XItvu31LbWSarkdx5gSIRSn6g4BmgNrVbUW\n0AawmbldEhcHTzwBL70EqakwbZrbiYKvfOny3NnoTkZ8PcLtKMaYAhSmeGSo6i4gTkTiVfVLnAv3\njIs6dXLGxXr0URg0KPamuO3btC+TV05m18FdbkcxxuSjMMVjr4iUB9KBySIyEufqcOOyCy90Lib8\n6iu45hqnPyRWVC1Xlc71O/PqYhuq3ZhIVJg+j3I4Q4HEATcDFYDJqro79PECF+t9Hvk5cgT69oX0\ndOeCwtoxMq/Sml1rSBmXgq+fj5NLnex2HGNiWtD7PFR1v6pmqWqGqo5X1ZGRWjhKqsREeOMN6N0b\nLrnEmS89FpxT+Rxa1mxpAyYaE4EK0/LYDxx9UiJQCtgfoms9iq0ktjxy+/JL6NrV6Qfp3dsZsTea\nLd6ymOvfv571fdZTKr6U23GMiVmhaHnkXCwIlAGuA14rRkYTQq1bOxcUvvWWUzyystxOVDxNqjeh\nbnJd3vn+HbejGGNyKdLYViKyXFUbhCBPsZX0lsdRv/0G110HSUnOuFgnneR2oqKb55tHz5k9WX3f\nauLjSsgY9caEWdBbHiLSKdftehF5BqcD3USwU06BTz5x5gNp2xZ+/dXtREXX6sxWVC1Xlff++57b\nUYwxfoU5VbcDcJX/dgXwO9AxlKFMcJQuDVOnQoMGzvhYW7a4nahoRITBKYMZmj6UbI2xC1qMiVKu\nDckeKnbY6q9U4fnn4bXXYPZsqFfP7USBU1WajGnCwJYDubbetW7HMSbmBHrYqsDiISKjci0qILnu\no6p9ixoylKx4FGziROjf3xkTq3lzt9ME7sPVHzJk/hC+7fWtTRZlTJAFs8/jO/+tNM4c4muBdUBD\nnFN2TZS59VYYNw6uvhpmznQ7TeCuPudqDmcdZs6GOW5HMabEK8x1Ht8ALf3ziiMipYCvVDXfKWLd\nZi2PEzs6P/qQIc78INFk2g/TGLV4FF/1+MpaH8YEUShG1a2IMyTJUeX960yUatIE5s+HoUOdAhJN\ntfb6+tez88BO5m2c53YUY0q0whSPZ4ClIjJBRCYAS4HhoY2VP3EMFZGRInKrGxliRZ06sHAhfPAB\n3Hdf9FxMGB8Xz8CUgQyZP8TtKMaUaIW5wnwc0Az4P+DfQDNVHR/iXAW5BqgOHMHmNy+2qlVh3jxY\nswZuuAEOHXI7UeHcfP7NrN+znkWbF7kdxZgSq8DiISL1/P82Ak4HNuH8wq4mIhcVZ6ci8raI7BCR\nlXnWp4rIahFZJyID8nlpXWCBqj4E3FOcDMZRoYJzMWFiIlxxRXRcTFgqvhQDWgxgaPpQt6MYU2Id\n71Td0araS0TS+HNgxByq2rrIOxVJwZkTZKKqnu9fFw+sAS4DtgDfAl1xJp66CHgeaA0cUdX3ReRd\nVe2Sz7atw7wIsrPhH/9wJpiaPRvOOMPtRMd3KPMQtUfWZlbXWTQ8vaHbcYyJekG7ziPURMQDzMxV\nPJoDT6hqqn/5EQBVfSbXa8oAo4CDOPOov57Pdq14FJEqvPACvPKKU0Dq13c70fG99PVLLNy8kPev\nf9/tKMZEvUCLR0IhNng9MEdV94nIYzjXeQxR1aXFyJmf6jiHxo7aDBxzOrCq/gH0PNGGvF4vHo8H\nj8eD1+vF6/UGNWisEoGHH3b6Qlq3di4mvOQSt1MV7M5Gd/LMgmf4ceeP1KsShZfNG+OitLQ00tLS\n8Pl8+Hy+gF9fmOs8Vqrq+SLSEhgCvAA8rqpNihI413Y9HNvy6ASkqmov/3I3oKmq9glwu9byCII5\nc+CWW2D8eLjySrfTFGxY+jBW71rNxGsnuh3FmKgWius8jp7EeRUwWlVn4UwIFWxbgBq5lmtgZ1S5\npm1bZ0rbHj3gvQgezPa+i+/jk3WfsGHPBrejGFOiFKZ4bBGRt4AuwMciclIhXxeoJUAdEfGISKJ/\nfzNCsB9TSM2aweefQ79+8HaEzgR7ykmncO/F9/LsgmfdjmJMiVKYInAD8Clwhar+CiQBDxdnpyIy\nFVgI1BWRTSLSQ1Uzgd7AHGAV8K6q/lic/Zjiu+ACSEuDJ5+Ef/3L7TT5u7/p/Xzw4wds+m3TiZ9s\njAmKQp35/eokAAAVf0lEQVRt5T+19mxVHSciVYDyqvpTyNMVgfV5hMYvv8Bllzn9IIMHR97c6P0/\n78+hzEOMbDfS7SjGRKWgn6orIv8EGgHnqGpdEakOvKeqLYqVNESseITOjh1w+eVOf8hzz0VWAdm+\nfzv1X63Pj/f9yGnlTnM7jjFRJxQd5tfizBx4AEBVt+AMjmhKmNNOcw5hpafD3XdH1nhYVctV5ebz\nb2bE1yPcjmJMiVCY4nFY9c+5P0WkbAjzmAhXqZLTib5unXMIKyPD7UR/erjFw4xZNobdB3e7HcWY\nmHfc4iHOhAmzRORNoKKI3AnMBcaEI5yJTOXLw8cfw7590KlT5AyoWPOUmlz3t+sY+Y31exgTasft\n8/AXj5XAA0Bb/+o5qvp5GLIVifV5hE9GhtP62LkTPvoIypVzOxFs2LOBZmObsaHvBiqUrnDiFxhj\ngCD3efh/C38H/KaqD/lvEVs4THiVKgWTJ8NZZzkd6Xv3up0IaleqTdvabXnt29fcjmJMTCvM2VZr\ngLOBjfg7zXHqygUhzlYk1vIIP1VnRN65c+Gzz5yOdTet2rmK1hNa81PfnyibaF10xhRGKE7V9eS3\nXlV9gQQLFyse7lCFp56CKVOcDvWaNd3N0/m9zrSs2ZJ+zfq5G8SYKBE1Q7KHihUPd730Erz8slNA\n6tRxL8eybcvoMLUDG/puoHRCafeCGBMlQnGdhzGF9sADzhXoXi98/717ORqe3pAGVRswfvl490IY\nE8Os5WFC4t13oW9fZ2Tepk1P/PxQ+HrT19z075tY23stpeJDMRC0MbHDWh4mInTp4ozEe9VV8MEH\n7mRoXqM5ZyWdxZSVU9wJYEwMs5aHCanvvoNrroGePeGxxyAuzH+ufPnzl9z98d2suncV8XHx4d25\nMVHEWh4mojRqBIsXO3Oid+kCBw6c+DXB5PV4qXxyZaavmh7eHRsT46x4mJA7/XRnQMWTT4aWLZ3h\n3cNFRBiUMohhXw3DWqTGBE9UFQ8RqSki/yciY0VkgNt5TOGddJIzH3q3bs4MhQsWhG/f7c5uR7zE\nM2vtrPDt1JgYF1XFAzgPmK6qdwAN3Q5jAiPiXIk+dixce234prYVEQamDGRI+hBrfRgTJK4UDxF5\nW0R2iMjKPOtTRWS1iKwroGXxDXCHiMzFmRrXRKF27WDePBg+3LkuJDMz9Pu8rt51/H74d+b+PDf0\nOzOmBHCr5TEOSM29QkTigVf86+sDXUWknojcIiIviUg14DbgCVVtA7QPc2YTRPXqOR3pP/wA7duH\nflDFOInj0ZaPMjR9aGh3ZEwJ4UrxUNV0IO+viybAelX1qWoGMA3oqKqTVPUBVd2K09roKyKvAz+H\nN7UJtqQk5yysv/3N6QdZsya0++t6flc2/rqRBb+EscPFmBiV4HaAXKoDm3ItbwaOuTZZVf8LXH+i\nDXm9XjweDx6PB6/Xi9frDWpQEzwJCc5YWGPGQEoKTJrkzJEekn3FJTCgxQCGpg/lk5s/Cc1OjIkS\naWlppKWl4fP58Pl8Ab/etYsE/aP1zlTV8/3LnYBUVe3lX+4GNFXVPgFu1y4SjFLp6XDDDfDww05f\niBT6cqXCO5x5mNojazOj6wwuOv2i4O/AmCgVzRcJbgFq5FqugdP6MCVESgosWgQTJ8Ltt8Phw8Hf\nR+mE0jx0yUMMSx8W/I0bU4JEUvFYAtQREY+IJAJdgBkuZzJhduaZ8NVXzvzorVvD9u3B30evi3qR\n/ks6q3auCv7GjSkh3DpVdyqwEKgrIptEpIeqZgK9gTnAKuBdVf3RjXzGXeXKwfvvwxVXQJMm8J//\nQFZW8LZfNrEs9ze9n+FfDQ/eRo0pYWxgRBPRPvjAmaFw+3bo0MG5uLBNG+eK9eL47dBv1B5Zm8W9\nFnNW0lnBCWtMFLOZBK14xKSffoKPPoIPP4Tly51WyTXXwJVXOqf8FsVjXzzG/w78jzc7vBncsMZE\nISseVjxi3s6dMHOmU0jS0pzJpq65Bjp2hDPOKPx2dh3cRd1Rdfn+nu85o0IALzQmBlnxsOJRouzf\nD5995hSSjz+G2rWdQnLNNc5V7Cc63fcfc/5Blmbxr9R/hSewMRHKiocVjxIrIwPmz3cKyYcfQpky\nfxaSswro1th+YCut3z+P9C6rqVzm1PAGNiZIypVzbsVhxcOKhwFUYelSp4jMmAE7dhT83H0p9yJH\nKlL+G7v2w0SnRx6Bfv2Ktw0rHlY8TIB8v/po9FYj1vdZT1KZIva+GxPlovkKc2Nc4anooUPdDryy\n+BW3oxgTNazlYQywZtcaUsal8NP9P1EusZgHj42JQtbyMKYIzql8Dl6PlzeX2DUfxhSGtTyM8Vux\nfQXtJrfjp/t/4qSEYl7CbkyUsZaHMUV0YdULaVStEW8vC9Pk6sZEMWt5GJPLos2LuHH6jazrs45S\n8aXcjmNM2FjLw5hiaHZGM2pXqs3klZPdjmJMRLPiYUweg1MGMyx9GFnZQRwH3pgYY8XDmDy8Hi+V\nT67M9FXT3Y5iTMSK2OIhIrVEZIyIvO9fLisiE0TkLRG5ye18JnaJCINSBjHsq2FY/5kx+YvY4qGq\nP6tqz1yrrgPeU9U7gatdimVKiCvrXEmcxDFr7Sy3oxgTkUJePETkbRHZISIr86xPFZHVIrJORAYU\nYlPVgU3++3Yw2oSUiDCw5UCGpA+x1ocx+QhHy2MckJp7hYjEA6/419cHuopIPRG5RUReEpFq+Wxn\nM1DDfz9iW0wmdlxX7zr2Hd7H3J/nuh3FmIgT8l/CqpoO7M2zugmwXlV9qpoBTAM6quokVX1AVbeK\nSCUReQNo6G+Z/BvoJCKvATNCnduY+Lh4Hm35KEPTh7odxZiIk+DSfnMfggKnVdE09xNUdQ9wd57X\n3V6YjXu9XjweDx6PB6/Xi9frLU5WU4J1Pa8rT6Q9wYJfFtCiZgu34xgTNGlpaaSlpeHz+fD5fAG/\nPixXmIuIB5ipquf7lzsBqaray7/cDWiqqn2CsC+7wtwE1RtL3mDGmhl8cvMnbkcxJmSi5QrzLfzZ\nf4H//maXshhzXLc1uI0VO1awdNtSt6MYEzHcKh5LgDoi4hGRRKAL1o9hItRJCSfx8CUPM2T+ELej\nGBMxwnGq7lRgIVBXRDaJSA9VzQR6A3OAVcC7qvpjqLMYU1R3NrqThZsWsnLHyhM/2ZgSwEbVNaaQ\nnlvwHMu2L2Nqp6luRzEm6KKlz8OYqHNP43uY+9Nc1uxa43YUY1xnxcOYQipfujx9mvRh2FfD3I5i\njOvssJUxAfj10K/UHlmbb3t9y1lJZ7kdx5igscNWxoRQxZMqck/je3jmq2fcjmKMq6zlYUyAdh3c\nRd1RdVl+93JqnlLT7TjGBIW1PIwJsconV6bnRT15fsHzbkcxxjXW8jCmCLbv3079V+vz33v/y+nl\nT3c7jjHFZi0PY8Kgarmq3HLBLbyw8AW3oxjjCmt5GFNEm/dt5oLXL2BN7zVUKVvF7TjGFIu1PIwJ\nkzMqnMEN597AS4tecjuKMWFnLQ9jisH3q49GbzViXZ91VCpTye04xhSZtTyMCSNPRQ8dz+nIqG9G\nuR3FmLCylocxxbRu9zouefsSNvTdQIXSFdyOY0yRWMvDmDCrk1yHK2pfwauLX3U7ijFhYy0PY4Jg\n1c5VtJ7Qmp/6/kTZxLJuxzEmYNbyMMYF9avUJ6VmCm9+96bbUYwJi4hueYhILWAQcIqqXi8iHYH2\nQAVgrKp+ns9rrOVhXLF8+3KunHwlG/puoEypMm7HMSYgMdXyUNWfVbVnruWPVPVO4G6cec+NiRgN\nqjagcbXGvL3sbbejGBNyYSkeIvK2iOwQkZV51qeKyGoRWSciAwLY5GDgleCmNKb4BrcazLMLnuVI\n1hG3oxgTUuFqeYwDUnOvEJF4nAKQCtQHuopIPRG5RUReEpFqeTcijmeB2aq6PBzBjQlEk+pNqFel\nHhOWT3A7ijEhFZbioarpwN48q5sA61XVp6oZwDSgo6pOUtUHVHWriFQSkTeABiLyCNAbaAN0FpG7\nwpHdmEA91uoxhn81nMzsTLejGBMyCS7uuzqwKdfyZqBp7ieo6h6c/o3cTngpr9frxePx4PF48Hq9\neL3e4mY1ptBa1mxJzVNqMmXlFG698Fa34xiTr7S0NNLS0vD5fPh8voBfH7azrUTEA8xU1fP9y52A\nVFXt5V/uBjRV1T7F3I+dbWVcN/enudz7yb2suncV8XHxbscx5oSi6WyrLUCNXMs1cFofxkS9v9f6\nO8llkpm+arrbUYwJCTeLxxKgjoh4RCQR59TbGS7mMSZoRITBrQYzJH0I2Zrtdhxjgi5cp+pOBRYC\ndUVkk4j0UNVMnA7wOcAq4F1V/TEceYwJh3Znt6N0fGk+Wv2R21GMCbqIvsK8KKzPw0SSD1d/yNPz\nn2ZJryWIFPpwsjFhF019HsbEvKvPuZojWUeYvX6221GMCSorHsaEUJzEMThlME/PfxprEZtYYsXD\nmBDrXL8ze//Yyxc/f+F2FGOCxoqHMSEWHxfPwJSBPDDnAZZtW+Z2HGOCws0rzI0pMbpd0I19h/dx\n5ZQr8Xq8POV9ijrJddyOZUyRWcvDmDCIkzh6N+nNuj7rOK/KeTQf25y7Z93N1t+3uh3NmCKx4mFM\nGJVLLMegVoNY03sNFUpX4LzXzmPA5wPY88cet6MZExArHsa4IPnkZJ67/DlW3rOSXw/9St1RdRmW\nPowDRw64Hc2YQrHiYYyLqleozpsd3mThHQtZsWMFZ486m1cXv2qTSZmIZ1eYGxNBlm5bysC5A1m7\ney1PtX6Krud1tVF5TVgEeoW5FQ9jIlCaL41H5z7KgSMHGPr3oVxV9yob3sSElBUPKx4mRqgqM9bM\nYNAXgzjlpFMY3mY4rc5s5XYsUwSqSmZ2Zr43RYmXeBLiEv5yi5O4sP3RYMXDioeJMVnZWUxZOYXH\n0x4nKzuLUvGliJd44iSO+Lh44iWe+Dj/sv/+8R4vU6oMZUuVpVxiOcqWKkvZxLJ/Ludzv2zisc+N\nE6erVFXJyM4gIyuDjOwMjmQdISPL/69//dH7Rx87ej+/X6JZ2VnHLmvWcR9XnJ9zQRCRY/4F/rLu\n6C/h3OsyszP/kvUvy3neX97HCyoKuW/Zmp1vcUiIS8j5P877fk/0uty3fzT/B3c3zjvpamCseFjx\nMDHqSNYRtuzbQpZmkZWdRbZm59zPUv+y/35Bj2dmZ3Io8xAHjhzgQMYB9h/Zf+z9jAMFP3bkAAcz\nDpIYn5jziy4hLoFScaVIjE+kVLz/37hSBd4/+ryjr0uISyA+zv9Xt/z5yzBn3dHlPH+ZHy2Gqoqi\nOf8Cf1l3vMfy5s8v+/EeT4hLyHk/+d2OPqeoLYij/6cnKk6VylQi+eTkYn2/rHhY8TAmZFSVQ5mH\niI+Lp1RcKeuHiSExMyS7iNQSkTEi8n6udWVF5FsRae9mNmNKKhGhTKkyJMYnWuEo4SK2eKjqz6ra\nM8/q/sC7buQpqdLS0tyOEDPsswwu+zzdFfLiISJvi8gOEVmZZ32qiKwWkXUiMqAQ27kcZ7ranaHK\nav7KfkCDxz7L4LLP013haHmMA1JzrxCReOAV//r6QFcRqScit4jISyJSLZ/tXAo0A24Ceom1mY0x\nxjUhH5JdVdNFxJNndRNgvar6AERkGtBRVZ8BJvnXVQKGAQ1EZICqDvav7w7stF5xY4xxT1jOtvIX\nj5mqer5/uTPQVlV7+Ze7AU1VtU8Q9mVFxRhjiiCQs63cmgwqZL/gA3nzxhhjisats622ADVyLdcA\nNruUxRhjTIDcKh5LgDoi4hGRRKALMMOlLMYYYwIUjlN1pwILgboisklEeqhqJtAbmINz+u27qvpj\nMfcT0Km/5vhExCci34vIMhFZ7HaeaJPfKeoiUklEPheRtSLymYhUdDNjNCng8/yniGz2f0eXiUjq\n8bZhHCJSQ0S+FJH/isgPItLXvz6g72dMDE/iP/V3DXAZziGxb4GuxS1IJZmI/Aw0UlWbH7UIRCQF\n2A9MzHWiyHPALlV9zv8HTpKqPuJmzmhRwOf5BPC7qo5wNVyUEZGqQFVVXS4i5YDvgGuAHgTw/YzY\nK8wDlHPqr6pmANOAji5nigV28kERqWo6sDfP6quBCf77E3B+YE0hFPB5gn1HA6aq21V1uf/+fuBH\noDoBfj9jpXhUBzblWt7sX2eKToH/iMgSEenldpgYcZqq7vDf3wGc5maYGNFHRFaIyFg7DBg4/2UU\nDYFvCPD7GSvFI/qPvUWeFqraEGgH3Oc/bGCCxH+Rq31vi+d1oBbQANgGvOhunOjiP2T1AXC/qv6e\n+7HCfD9jpXjYqb9Bpqrb/P/uBP4P59CgKZ4d/uPNiMjpwP9czhPVVPV/6geMwb6jhSYipXAKxyRV\n/dC/OqDvZ6wUDzv1N4hE5GQRKe+/Xxa4Alh5/FeZQpgBdPff7w58eJznmhPw/4I76lrsO1oo/nEB\nxwKrVPVfuR4K6PsZE2dbAYhIO+BfQDwwVlWHuxwpaolILZzWBjijEEy2zzMw/lPULwUq4xw/fhz4\nCHgPqAn4gBtU9Ve3MkaTfD7PJwAvziErBX4G7sp1zN4UQERaAvOB7/nz0NSjwGIC+H7GTPEwxhgT\nPrFy2MoYY0wYWfEwxhgTMCsexhhjAmbFwxhjTMCseBhjjAmYFQ9jjDEBs+JhjDEmYFY8jDHGBMyK\nhzEBEJGyIvKxiCwXkZUi0l9EPvA/1lFEDopIgoicJCIb/Otri8hs/wjF80XkHP/6KiIyXUQW+2+X\n+Nf/U0QmichC/8Q8Pd17x8bkL8HtAMZEmVRgi6q2BxCRCsBd/sdScMZXagKUAhb517+FM3TGehFp\nCrwGtAFeBl5S1QUiUhP4FKjvf815QDOgHLBMRD4+OlilMZHAiocxgfkeeEFEngFmqepXIrJBRP4G\nXAyMAFrhjLGW7h9Y8hLgfWc8OgAS/f9eBtTLtb68//kKfKSqh4HDIvIlTkH6KPRvz5jCseJhTABU\ndZ2INATaA0NEZC7OIHNXAhnAXJxZ2OKAh3CKyF7/3Ch5CdBUVY8cs1LynRwvO2hvwpggsD4PYwLg\nHwb8kKpOBl4ALgLSgX7AQlXdBSQDdVX1v6q6D/hZRDr7Xy8icoF/c58BfXNtu8HRu0BHESktIsk4\no8d+G/p3Z0zhWcvDmMCcDzwvItk4LY27ceaAPhWnBQKwgmOn8LwZeF1EBuP0hUzFOfzVF3hVRFbg\n/CzOA+7FOWz1PfAlzhDkT6nq9hC/L2MCYkOyGxNhROQJYL+q2rSqJmLZYStjIpP9VWcimrU8jDHG\nBMxaHsYYYwJmxcMYY0zArHgYY4wJmBUPY4wxAbPiYYwxJmBWPIwxxgTs/wFmo7yhh1doXgAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110fd4e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(fixedPoint[...,0], fixedPoint[..., 1], label=\"fixed point\")\n",
    "plt.semilogy(newton[...,0], newton[..., 1], label=\"newton\")\n",
    "plt.ylabel(\"residual\")\n",
    "plt.xlabel(\"sweep\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## coupled"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For \n",
    "$$ \\mathbf{F}(\\mathbf{x}) \\equiv \\left\\{\\begin{aligned}\n",
    "    \\frac{\\partial A}{\\partial t} &= \\nabla\\cdot\\left[\\left(3 + A^3\\right)\\left(2 + B^2\\right)\\nabla A\\right] \\\\\n",
    "\\frac{\\partial B}{\\partial t} &= \\nabla\\cdot\\left[A^5\\nabla B\\right]\n",
    "\\end{aligned}\\right.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A = fp.CellVariable(mesh=mesh, name=\"A\", hasOld=True)\n",
    "B = fp.CellVariable(mesh=mesh, name=\"B\", hasOld=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aface = A.faceValue\n",
    "Bface = B.faceValue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Aeq = fp.TransientTerm(var=A) == fp.DiffusionTerm(coeff=(3 + Aface**3)*(2 + Bface**2), var=A)\n",
    "Beq = fp.TransientTerm(var=B) == fp.DiffusionTerm(coeff=A**5, var=B)\n",
    "eq = Aeq & Beq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A.constrain(1., where=mesh.facesLeft)\n",
    "B.constrain(2., where=mesh.facesRight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A.value = 0.\n",
    "B.value = 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fixedPoint = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A.updateOld()\n",
    "B.updateOld()\n",
    "for sweep in xrange(20):\n",
    "    res = eq.sweep(dt=1000.)\n",
    "    fixedPoint.append([sweep, res])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fixedPoint = fp.numerix.array(fixedPoint)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAEACAYAAABPpeiSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHnpJREFUeJzt3Xl0VeW9//H3NwNhDDFgAhJKFAEFB8ApTiXiRB1bpVWq\nt7Tan14VbbW3k7eteO+qq+2yt7WD1k5cUJy9i6JSkKpR0IpIARECJkI0gRAiCVOAkOH5/fEcJIac\nJOQM+5zk81prrzPt7P015nx4nmc/e29zziEi0paUoAsQkcSlgBCRsBQQIhKWAkJEwlJAiEhYCggR\nCavdgDCz4Wb2mpmtNbP3zeyuMOv9xsxKzGy1mU2ITakiEm9pHXzeANztnFtlZv2BFWa22DlXfHAF\nM7sMON45N8rMzgIeAQpiV7KIxEu7LQjn3Fbn3KrQ8z1AMXBMq9WuAmaH1lkGZJlZbgxqFZE46/QY\nhJnlAxOAZa0+GgaUt3hdAeRFWpiIBK9TARHqXjwHfCvUkjhslVavNX9bpBvoaAwCM0sHngced87N\na2OVzcDwFq/zQu+13o5CQyQgzrnW/4h3SkdHMQz4C7DOOffrMKvNB74WWr8A2OGcqwpTZNIs9913\nX+A1qObEW5KtXuci+3e5oxbEucCNwHtmtjL03r3A50Jf+EedcwvM7DIzKwXqgG9EVJGIJIx2A8I5\nt5ROjFM452ZErSIRSRiaSRlGYWFh0CUcMdUce8lWb6Qs0j5Kp3dk5uK1LxE5xMxwXRyk7PAoRqz5\ncdDuQyEo3UngAQHd50vV3cJORGMQIhKWAkJEwlJAiEhYCggRCUsB0QmFhYVkZ2dz4MCBoEsRiSsF\nRAfKyspYsmQJKSkpzJ8/P+hyROJKAdGBOXPmcPbZZzN9+nRmz54ddDkinbaodBErK1d2vGI7FBAd\nmDNnDjfeeCM33HADixYtYtu2bUGXJNIpc9fM5b2q9yLaRsIHhFl0lq5YunQpH3/8MV/5yleYOHEi\nI0eO5Iknnojuf6BIjFTvrebofkdHtI2EDwjnorN0xezZs7nkkkvIzs4GYNq0aepmSNLYVreNo/tG\nFhAJMdU6Ee3bt49nnnmG5uZmhg4dCkB9fT07duzgvffe45RTTgm4QpH2Vdf1gBZEUObNm0daWhrF\nxcWsXr2a1atXU1xczPnnn8+cOXOCLk+kXc4538WIsAWhgAhjzpw53HTTTeTl5ZGTk0NOTg65ubnM\nmDGDJ554gubm5qBLFAmrrqGOFEuhX69+EW0n8OtBhM5Vj0sNsdad/lskuW2s3cjk2ZMp+3ZZRNeD\nUAtCpBuKxvgDKCBEuqXqvdXk9MuJeDsKCJFuqLou8gFKUECIdEvRmAMBCgiRbikasyhBASHSLWkM\nQkTC0hiEiIS1rW6buhgi0rZoTLMGBURY+fn59O3blwEDBpCdnc0VV1xBRUVF0GWJdEp1ncYgYsrM\nePHFF9m9ezeVlZXk5uZy5513Bl2WSIfqDtQBRHweBiggOiUjI4Nrr72WdevWBV2KSIeiNf4ACoh2\nHTzxau/evTz99NOcffbZAVck0rFojT9AElwwxu6Pzv0u3X1Hdpalc44vfvGLpKWlUVdXR05ODgsX\nLoxKLSKxFK0TtSAJAuJIv9jRYmb87W9/Y/LkyTjnmDdvHpMmTWLdunXk5uYGUpNIZ0RrkhSoi9Ep\nZsaXvvQlUlNTefPNN4MuR6Rd0ToPA5KgBRGkg2MQzjnmz59PbW0tJ554YsBVibQvWrMoQQHRriuv\nvJLU1FTMjPz8fObMmaOAkIRXvbeaE4+Ozt+pAiKMTZs2BV2CSJdoDEJEwormGIQCQqSbqdpTRW7/\n6BxpU0CIdCPOOarqqhjSf0hUtqeAEOlGavfX0je9L73TekdlewoIkW6kcndl1FoPoIAQ6Va27tnK\n0P5Do7a9hDjMaRad8y1Eerqte7ZGtQUReEDoVnUi0VO5R10MEQkj2l0MBYRINxLtLkaHAWFmfzWz\nKjNbE+bzQjPbaWYrQ8uPoladiByRyj2VDB0Q30HKWcBvgTntrPO6c+6q6JQkIl0V9xaEc24JUNvB\najoMIZIAEnEehAPOMbPVZrbAzMZGYZsicoTqG+vZc2AP2X2yo7bNaBzm/Bcw3Dm318y+AMwDRre1\n4syZMz99XlhYSGFhYRR2LyIAVXX+JK03Xn+DoqKiqGzTOjMPwczygReccyd3Yt1NwGnOuZpW7zvN\neRCJnWUVy5jx9xks/3/LP/O+meGc69IwQMRdDDPLtdBUSDM7Ex86NR38mIhEWbTnQEAnuhhm9iQw\nCRhsZuXAfUA6gHPuUWAqcJuZNQJ7geujWqGIdEq0j2BAJwLCOTetg89/D/w+ahWJSJdU7qmMegtC\nMylFuolYtCAUECLdhAJCRMKK9jRrUECIdBtqQYhIm5xzCggRaVvNvpqoXqz2IAWESDewefdm8jLz\nor5dBYRIN1Cxq0IBISJtq9hVQd4ABYSItKF8Z7laECLStord6mKISBgagxCRsBQQItIm55zGIESk\nbbvqd2FmZGZkRn3bCgiRJHewexGLe9wqIESSXKzGH0ABIZL0FBAiElasZlGCAkIk6akFISJhxWoW\nJSggRJKeWhAiEpYCQkTatOfAHvY37o/qDXtbUkCIJLHNuzbHbJIUKCBEklrFrgqGDRgWs+0rIESS\n2Ec7P2JE1oiYbV8BIZLEynaUkT8wP2bbV0CIJLGyHWXkZ+XHbPsKCJEkVrajjGOPOjZm21dAiCSx\nTTs2qQUhIoc70HSAqj1VMZskBQoIkaRVvrOcYwYcQ1pKWsz2oYAQSVKxHqAEBYRI0lJAiEhYZTvK\nODYrdkcwQAEhkrTKdqoFISJhqIshImFtqo3tHAhQQIgkpfrGerbVbWNYZuzO5AQFhEhSKt9VzrDM\nYTGdAwEKCJGkFI/xB1BAiCSleBziBAWESFKKxwAlKCBEklJpbSnHZx8f8/0oIESSUMn2EkZlj4r5\nfhQQIknGOUdJTUlitCDM7K9mVmVma9pZ5zdmVmJmq81sQnRLFJGWquqqyEjN4Kg+R8V8X51pQcwC\npoT70MwuA453zo0CbgEeiVJtItKGku0ljBoU++4FdCIgnHNLgNp2VrkKmB1adxmQZWa50SlPRFor\nqYnP+ANEZwxiGFDe4nUFELtrYIn0cPEaoASI1jzN1vf9cm2tNGHCTAoKIDcXCgsLKSwsjNLuRXqO\nkpoSpo6dGvbzoqIiioqKorKvaATEZmB4i9d5ofcOM3XqTB5+GEaPhvHjoakJUlOjUIFID9JRF6P1\nP773339/l/cVjS7GfOBrAGZWAOxwzlW1teJ//ids2gS33AI//zkcfzw8+CDUtjfCISKfcs5RWlOa\nOIOUZvYk8BYwxszKzewmM7vVzG4FcM4tADaaWSnwKHB7e9vr1QumTYN//hOefhpWr4bjjoNbb4U1\nYQ+kighA5Z5K+vfqT2ZGZlz212EXwzk3rRPrzOjKzs88Ex57DKqq4I9/hClTYNQomDEDrr4a0tO7\nslWR7iueA5SQIDMpc3Phxz+GsjK4/XZ46CE49lj4r/+CysqgqxNJHCU18ZsDAQkSEAelp8NXvgJL\nlsCCBbBlC4wdC9ddB6+/Dq7NYyMiPUePbEG05ZRT4A9/8K2K88+H226DcePgN7+BHTuCrk4kGPGc\nJAUJHBAHDRzoxyTWroVHHoG33vLdj5tvhuXL1aqQnmX9J+sZPWh03PaX8AFxkBlMmgRPPQUbNvjB\nzOuug9NOg0cfhV27gq5QJLYONB1g045NjBk8Jm77TJqAaCknB37wAygthZ/9DF5+GUaMgG9+E5Yt\nU6tCuqfSmlI+N/Bz9E7rHbd9JmVAHJSSApdcAs8/D8XFfuLVDTf4WZq//S3U1ARdoUj0rN22lrFH\nj43rPpM6IFoaMsS3Kj74AH71K3j7bT8B66tfhVdegebmoCsUicy66nWMHayAiEhKCkyeDHPnwsaN\ncPbZ8J3vwMiRcP/98NFHQVco0jVrq9cyLmdcXPfZ7QKipexsuPNOWLnSd0Oqq2HiRLjoInj8cdi7\nN+gKRTpvXfW6uHcxzMVpRM/MXLz21Z79+2H+fJg1y3dDrr0Wpk+H887zR0pEElFDUwOZP8uk5ns1\n9Envc0Q/a2Y457r0192tWxBt6d3bz9b8+9/93IrRo+Hf/913QWbOhA8/DLpCkcOV1pSSl5l3xOEQ\nqR4XEC0dcwx873vw/vvwzDP+tPNzzvHLI4/A9u1BVyjiBdG9gB4eEAeZwemn+5PEKirg3nvhjTf8\nUZCrr/bhsW9f0FVKT7a2ei3jjo7vACUoIA6Tng5XXAFPPgnl5XDNNfCnP/nWxvTpsGgRNDYGXaX0\nNGpBJKDMTB8KixfDunUwYQL85CcwbBjccQcsXar5FRIfQbUgetxRjGj48EN/TsjTT/txiy9/2Q98\nnnWWjoRI9B08grH9e9vpm973iH9eRzHibORIf33N997zXY4BA+DrX/dnmf7Hf8A77+h8EIme4k+K\nyc/K71I4REoBEaGxY/0MzeJiP7+iTx/42tcgPx/uucefnq5uiERiZeVKJgwJ5o6WCogoMfMXufnv\n//Zh8dJLvmVxyy2Ql+fHLF55BRoagq5Uks3KrSsZP2R8IPtWQMSAGZx0km9ZvP8+vPYaDB/uTyYb\nMsQPfM6bp6ne0jmrtq4KrAWhQco4Ky/34TBvnr8i1gUX+LkWV1zhr3Mh0pJzjqN+fhSld5UyuO/g\nLm0jkkFKBUSAamr8xXnnz/cXvRk7Fq680i/jxumIiMDG2o1M+t9JlN9d3vHKYUQSENG6N6d0QXY2\n3HijX+rr/ZW758/3rQkzuPxyv1xwgT+HRHqeILsXoDGIhJGR4a+O9bvf+dsTzp/vBzcfeMB3Pa68\nEh5+2H8mPcfKyuAGKEEBkZDM4OST/aDmkiX+0v833uivt1lQACecAN/+tj8jVQOd3dvKrcEd4gSN\nQSSd5mZ/AZxFi2DhQv+8oAAuvtgvp57qr6ol3UPe/+Sx9Kal5Gfld3kbGqTswXbt8odRFy/2S20t\nXHihv2rWhRf6CVuSnKrrqhn9u9HUfK8Gi2DEWoOUPVhmpj9MevXV/vXHH/sJWYsX+9PW+/f3g5yT\nJ/vHoUODrVc6b0XlCiYOnRhROERKLYhuzDl/1azXXoNXX/VHSXJyoLDQ34Ro0iR/GrskpplFM2lo\nauCnF/40ou2oiyGd0tTkTzArKvJhsWSJP9T6+c/7+5+ed54/EU3zLxLDlMencMcZd3DlmCsj2o4C\nQrqkudm3MN54w4fFkiX+vfPOg3PP9cv48f4iOhJfza6ZQb8YxIYZG8jpF9kUWwWERIVz/pDq0qX+\nLNQ33/T3Fpk40V+ns6DAL0OGBF1p97f+k/Vc/sTlfHhX5FdRVkBIzOzc6a9v8dZbfh7G22/7s1TP\nOssvZ57pA6Rfv6Ar7V5mrZzFPzb9g7nXzI14WzqKITEzcOChORbgWxklJT4s3nkHnn0W1qzxF/g9\n4wx/t/XTT/envveJ7xXau5W3K96mYFhB0GWoBSGRO3DAh8Ty5bBiBbz7LmzY4G+mPGHCoeXUUyEr\nK+hqk8OpfziVP1/5Z84YdkbE21IXQxJOfb0PjVWr/GzPlSv968GDfevi1FP948kn+yBJTQ264sSx\nu343Q345hNrv19IrtVfE21MXQxJORobvapx++qH3mpv9BX9Xr/bL44/7w65bt/rzS046yS/jxvlT\n30eM6JnTxt/d8i7jh4yPSjhESgEhcZOSAqNG+WXq1EPv79njD7e+/75/fOUVf5uBmhp/a8QTT/TL\nCSfAmDG+xdE3/tdvjZu3yt9KiPEHUBdDEtiuXbB+vb/G54YN/vGDD/yh15ycQ2EzapQPjZEj/WBp\nsg+OXvzYxdx15l0RT5A6SGMQ0qM0NcFHH/mjKSUlUFrqHz/80M/jGDTI34Kg5ZKf75e8vMSe+HWg\n6QCDfjGI8rvLyeodnRFdjUFIj5Ka6lsKxx0Hl1762c+ammDLFt/K2LjRB8Zrr/lAKSuDykrf+vjc\n5/yFhIcP96GRl+fvmDZsmJ8I1iug7v87m99hzKAxUQuHSCkgpFtJTT30xZ806fDPGxth82Z/8eCD\ny6ZNfrr5li3+s6oqP/9j6FC/DBkCubmHlpwcOProQ0tGRvTqLyorojC/MHobjJACQnqUtDR/dGTE\niPDrNDfDtm3+6EplpX+sqvJ3fv/Xv/zz6mq/fPKJD4hBg/ySnQ1HHeUfs7L884ED/ZKZeWjp3//Q\n0qfPoaM1r5W9xj0F98Tnl9EJGoMQiYBzfjB1+3Z/1GX7dn/RnpoaP029ttY/7tgBu3f7dXft8kdu\ndu+GujrYv9+HRJ/+9dR8czDHzqugX+pAMjL4dOnVy4+dpKf7kEtL862llJRDj2afXQCmTYNJkzQG\nIRIIs0MthOOO69o2mpth3z54pXQZP156As/830D27/eTzfbv93djq6/3jw0NfpylocH/XFOTX5zz\nr507dF9Y53xLJhIKCJGApaT4k91W7Sji0tEXMGZM0BUd0uE8NTObYmbrzazEzL7fxueFZrbTzFaG\nlh/FplSR7u0fG//BBfkXBF3GZ7TbgjCzVOB3wEXAZmC5mc13zhW3WvV159xVMapRpNur3VfLyq0r\nE+oIBnTcgjgTKHXOlTnnGoCngKvbWE8XKROJwOKNi/n8iM/TJz2xpoF2FBDDgJY3BawIvdeSA84x\ns9VmtsDMxkazQJGeYEHJAi4fdXnQZRymo0HKzhyX/Bcw3Dm318y+AMwDRre14syZMz99XlhYSGFh\nYeeqFOnGml0zfy/9O/dNui8q2ysqKqKoqCgq22p3HoSZFQAznXNTQq9/CDQ7537ezs9sAk5zztW0\nel/zIETasHzzcqbPm866O9bFZPuRnIvRURfjXWCUmeWbWS/gOmB+q53nWujOHmZ2Jj50ag7flIi0\nZUHJAi4bdVnQZbSp3YBwzjUCM4BFwDrgaedcsZndama3hlabCqwxs1XAr4HrY1mwSHezoDQxxx9A\nU61FAlW+s5zxj46n8juVMbuCVCy7GCISQ88XP8/VY65OiMvLtUUBIRKgZ9c9y5fHfjnoMsJSQIgE\npHxnOes/Wc+Fx10YdClhKSBEApLo3QtQQIgE5pm1zyR09wIUECKB+Hjnx2zYviGhuxeggBAJxOxV\ns7l+3PUJ3b0ABYRI3DW7ZmatmsVNE24KupQOKSBE4uz1stcZkDGAiUMnBl1KhxQQInH211V/5eYJ\nN2OW+JdRUUCIxNGO/Tt4YcML3HDyDUGX0ikKCJE4mvveXC4ZeQmD+g4KupROUUCIxEmza+ahZQ9x\n11l3BV1KpykgROLkpQ9eYmDvgZw7/NygS+k0BYRInPzq7V9xd8HdSTE4eZACQiQOVm1dRUlNScJP\nrW5NASESBw++9SAzzphBemp60KUcEQWESIwVVxfz8ocvc9sZtwVdyhFTQIjE2P2v3889Z99DZkZm\n0KUcMQWESAytqVpDUVkRM86cEXQpXaKAEImhma/P5LvnfJf+vfoHXUqXKCBEYuT1std5d8u7STn2\ncJACQiQGGpsbuWvhXTx48YP0Te8bdDldpoAQiYE/rvgj2X2ymTp2atClRKSjm/eKyBHaVreNmUUz\neeVrryTVrMm2qAUhEkXOOW5/6Xa+Mf4bnJx7ctDlREwtCJEoenrt0xR/Uszj1zwedClRoYAQiZLK\n3ZV8a+G3eHHai/RO6x10OVGhLoZIFDQ2NzLt+WncfvrtnDHsjKDLiRoFhEgU3PvKvfRO682PPv+j\noEuJKnUxRCL03LrneGbtM6y4ZQWpKalBlxNVCgiRCLz58Zvc9tJtLLxhYdJcZ/JIqIsh0kXF1cVc\n88w1PPalxzjtmNOCLicmFBAiXVCyvYRLH7+UX1z0C6YcPyXocmJGASFyhIqri7lg9gXcN+k+po+f\nHnQ5MaWAEDkC72x+hwvnXMgDFz7AzRNvDrqcmNMgpUgnPbv2WW5fcDt/ueovXDXmqqDLiQsFhEgH\nGpoa+MlrP2Humrks/rfFjB8yPuiS4kYBIdKOsh1lfPX5r5LVO4sVt6zg6H5HB11SXGkMQqQNTc1N\nPPT2Q5z+x9O55sRrePGrL/a4cAC1IEQOs/Tjpdy96G76pPXhzZveZMzgMUGXFBgFhEjIe1XvMbNo\nJisqV/DA5AeYdvI0UqxnN7IVENKjOedY+vFSfvnPX7Js8zK+c/Z3mHvNXPqk9wm6tISggJAeqXZf\nLU+9/xSPrniUfY37uPPMO3ny2icVDK2Ycy4+OzJz8dqXSFt21e/ixQ9e5Ll1z/Hqple59PhLuWn8\nTVw88uJu3ZUwM5xzXbo4pgJCuq3G5kZWbV3Fq5teZWHpQpZvWc6kEZO49sRr+eIJX+SoPkcFXWJc\nxDQgzGwK8GsgFfizc+7nbazzG+ALwF7g6865lW2so4CQmHHOsXn3ZlZWrmT5luUs27yMZRXLyMvM\nozC/kCnHT6EwvzBp73AViZgFhJmlAhuAi4DNwHJgmnOuuMU6lwEznHOXmdlZwEPOuYI2tpVUAVFU\nVERhYWHQZRyRnlDz3oa9fLTjIz6s/ZCS7SV8sP0D1lavZW31WtJS0pgwZAKnDT2NgrwCCvIKoj53\nIRl/x5EEREeDlGcCpc65stCOngKuBopbrHMVMBvAObfMzLLMLNc5V9WVghJFMv4hJGvN555/LrX7\na6nZV8Mnez+huq6abXXbqKqronJ3JZV7KqnYVUH5rnJ27t/JiKwRHJt1LKOyRzH26LFMHTuVcTnj\nyO2XG/P7UCTj7zgSHQXEMKC8xesK4KxOrJMHJHVA9CTOOZpcE03NTTQ2N9LY3EiT888bmhpoaG74\nzOOBpgOfLvVN9dQ31lPfVM/+xv2fLnsb9rKvYR91DXXsbdhLXUMdew7sYXf9bnYf2M2u+l3s2L+D\n6iXV/PSBn5LVO4vsPtkM7juYwX0Hk9svl5x+OZyUcxIXj7yYvMw8hmcOJ7d/brceUEw0HQVEZ/sE\nrWO7zZ+7aM5Fndxcyw11vVvSXpemve065yhbVUbR/xYdtm7rbR787OD7ba3bch2H+8w2Dr5u+dh6\n3daPza75sOfNrpnat2uZ9etZNLvmw5am5qZDz13Tp+8dfJ5iKaRaKump6aRaKmkpaaSnpvvHlPRP\nX6enpJORlkF6Sjq9UnuRkZZBRmoGvdN6k5GWQZ+0PvRO603f9L70SetDbr9c+vXqR9/0vgzoNYD+\nvfqTmZHJgIwBZPXO4vc1v+eBHz2Q9Heg6q46GoMoAGY656aEXv8QaG45UGlmfwCKnHNPhV6vBya1\n7mKYWfIMQIh0M7Eag3gXGGVm+cAW4DpgWqt15gMzgKdCgbKjrfGHrhYoIsFpNyCcc41mNgNYhD/M\n+RfnXLGZ3Rr6/FHn3AIzu8zMSoE64Bsxr1pE4iJuE6VEJPnEfDjYzKaY2XozKzGz78d6f11hZsPN\n7DUzW2tm75vZXaH3s81ssZl9YGYvm1lW0LW2ZGapZrbSzF4IvU70erPM7DkzKzazdWZ2VhLUfHfo\nb2KNmT1hZhmJVLOZ/dXMqsxsTYv3wtZnZj8MfRfXm9klHW0/pgERmmj1O2AKMBaYZmYnxnKfXdQA\n3O2cGwcUAHeE6vwBsNg5Nxp4JfQ6kXwLWMeho0aJXu9DwALn3InAKcB6ErhmMxsG3Amc5pw7Gd/N\nvp7EqnkW/vvVUpv1mdlY/Dji2NDPPGzWwTFj51zMFuBsYGGL1z8AfhDLfUap7nn42aPrgdzQe0OA\n9UHX1qLGPOAfwAXAC6H3ErnegcDGNt5P5JqHAR8DR+HH614ALk60moF8YE1Hv1Pgh8D3W6y3ECho\nb9ux7mK0NYlqWIz3GZHQEZsJwDL8L/ngEZkqIDegstryK+C7QHOL9xK53mOBajObZWb/MrM/mVk/\nErhm59xm4Jf4kNiCP0K3mASuOSRcfcfgv4MHdfh9jHVAJNUIqJn1B54HvuWc293yM+cjNyH+e8zs\nCmCb8yfFtXn4OJHqDUkDJgIPO+cm4o94faZpnmg1m9lR+FMJ8vFfrv5mdmPLdRKt5tY6UV+7tcc6\nIDYDw1u8Hs5nEyxhmFk6Phwec87NC71dZWZDQp8PBbYFVV8r5wBXmdkm4Elgspk9RuLWC/7/e4Vz\nbnno9XP4wNiawDVfBGxyzm13zjUC/4fvNidyzRD+76D19zEv9F5YsQ6ITydamVkv/ADJ/Bjv84iZ\nn+f7F2Cdc+7XLT6aDxy8t9p0/NhE4Jxz9zrnhjvnjsUPmr3qnPs3ErReAOfcVqDczEaH3roIWIvv\n1ydkzcBHQIGZ9Qn9jVyEHxRO5Joh/N/BfOB6M+tlZscCo4B32t1SHAZQvoA/ZbwU+GGQgznt1Hge\nvi+/ClgZWqYA2fiBwA+Al4GsoGtto/ZJwPzQ84SuFzgVf8mA1fh/jQcmQc0z8Wcvr8GftZyeSDXj\nW5BbgAP48b5vtFcfcG/ou7geuLSj7WuilIiEpfNmRSQsBYSIhKWAEJGwFBAiEpYCQkTCUkCISFgK\nCBEJSwEhImH9f3Qr+BbUvej5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110fe1510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "viewer = fp.Viewer(vars=(A, B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Jacobian is then\n",
    "$$\\delta \\mathbf{F} \\equiv \\left\\{\\begin{aligned}\n",
    "\\frac{\\partial\\, \\delta A}{\\partial t} \n",
    "&= \\nabla\\cdot\\left[\\delta A \\, 3 A^2\\left(2 + B^2\\right)\\nabla A\\right] \n",
    "+ \\nabla\\cdot\\left[\\left(3 + A^3\\right)\\left(2 + B^2\\right)\\nabla \\delta A\\right] \n",
    "+ \\nabla\\cdot\\left[\\delta B\\, 2 B \\left(3 + A^3\\right)\\nabla A\\right] \\\\\n",
    "\\frac{\\partial \\, \\delta B}{\\partial t} \n",
    "&= \\nabla\\cdot\\left[\\delta A\\,5 A^4\\nabla B\\right]\n",
    "+ \\nabla\\cdot\\left[A^5\\nabla \\delta B\\right]\n",
    "\\end{aligned}\\right.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dA = fp.CellVariable(mesh=mesh, name=r\"$\\delta A$\", hasOld=True)\n",
    "dB = fp.CellVariable(mesh=mesh, name=r\"$\\delta B$\", hasOld=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dAeq = ((fp.TransientTerm(var=dA) \n",
    "         == fp.ConvectionTerm(coeff=3*Aface**2*(2+Bface**2)*A.faceGrad, var=dA)\n",
    "         + fp.DiffusionTerm(coeff=(3 + Aface**3)*(2 + Bface**2), var=dA)\n",
    "         + fp.ConvectionTerm(coeff=2*Bface*(3 + Aface**3)*A.faceGrad, var=dB))\n",
    "        + fp.ResidualTerm(equation=Aeq))\n",
    "dBeq = ((fp.TransientTerm(var=dB)\n",
    "         == fp.ConvectionTerm(coeff=5*Aface**4*B.faceGrad, var=dA)\n",
    "         + fp.DiffusionTerm(coeff=Aface**5, var=dB))\n",
    "        + fp.ResidualTerm(equation=Beq))\n",
    "deq = dAeq & dBeq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dA.constrain(0., where=mesh.facesLeft)\n",
    "dB.constrain(0., where=mesh.facesRight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A.value = 0.\n",
    "B.value = 0.\n",
    "dA.value = 0.\n",
    "dB.value = 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "newton = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A.updateOld()\n",
    "B.updateOld()\n",
    "dA.updateOld()\n",
    "dB.updateOld()\n",
    "for sweep in xrange(20):\n",
    "    res = deq.sweep(dt=1000.)\n",
    "    A.value = A.value + dA.value\n",
    "    B.value = B.value + dB.value\n",
    "    newton.append([sweep, res, max(abs(dA)), max(abs(dB))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "newton = fp.numerix.array(newton)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAD7CAYAAACWhwr8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHdJJREFUeJzt3Xl4VdW9//H3NwORkDAETJhlUsQqVqQIDpeoqIBD7Q9/\nWofW6m3rrdX29rn2V/X2Wnyu3sfrvXWo2trbQcGKgNUi9KoUhShYwYlJIUAYZAgmhiFkgAxk/f5Y\niUY8h4Rkn7PPST6v59nPyRmy9xfkfFxr7bXXNuccIiKRpIRdgIgkLgWEiESlgBCRqBQQIhKVAkJE\nolJAiEhUafE6kJnpfKpISJxz1pbfi2sLwjmXNNsvfvGL0GtQzYm3JVO95z99Pq9tfq1d31l1MUQ6\nqLLqMvpk9mnXPhQQIh2UAiKG8vPzwy7hmKnm2EuWep1z7Dm4h96Zvdu1H3MuPmOHZubidSyRzq6i\npoJ+v+xH5d2VmBmujYOUcTuLEY1Zm+pOWApBSQRBdC8gAQICOs6XqqOFnSSvoAJCYxAiHVAQ4w+g\ngBDpkNSCEJGoyqrL6NNVASEiEagFEUf5+fnk5ORQW1sbdikirbKnWmMQcbFt2zaWLl1KSkoK8+fP\nD7sckVYpO6gWRFzMnDmTCRMmcOONNzJjxoywyxFplQ41DyKRzZw5kzvuuINx48Yxfvx4SktLyc3N\nDbsskaPqNGMQZsFsbbFs2TK2b9/O1VdfzZgxYxg+fDizZs0K9g8oEgN7qvfQu2snGINwLpitLWbM\nmMHFF19MTk4OANdee626GZLwnHOUVZcFMkipLkYUBw8eZO7cuTQ0NNCvXz8Aampq2L9/P2vWrGH0\n6NEhVygSWUVtBRlpGRyXdly795XwLYiwzJs3j7S0NNavX8/q1atZvXo169ev57zzzmPmzJlhlycS\nVVDjD6CAiGrmzJncfPPNDBw4kNzcXHJzc8nLy+O2225j1qxZNDQ0hF2iSERBjT9AAqwH0Xitelxq\niLWO9GeR5PXKpld4dMWjvHrDqwDtWg9CLQiRDkZdDBGJSgEhIlHtORjcGIQCQqSDUQtCRKJSQIhI\nVAoIEYkqqGnWoIAQ6XBKq0rJ65YXyL4UECIdSH1DPfsO7VMLItaGDBlCZmYm2dnZ5OTkcNlll7Fz\n586wyxI5qrLqMnod14u0lGCuw1RARGFm/PWvf6WiooLdu3eTl5fH7bffHnZZIkdVUllCXlYw3QtQ\nQLRKRkYG06ZNY926dWGXInJUQY4/gALiqJouvKqurmbOnDlMmDAh5IpEjq6kqoTcbsEtiZjwC8bY\nvcHc79L94tiusnTOceWVV5KWlkZVVRW5ubm8+uqrgdQiEisllSWBtiASPiCO9YsdFDPjpZde4oIL\nLsA5x7x585g4cSLr1q0jLy+4/wAiQSqtKtUYRLyZGd/4xjdITU3lrbfeCrsckag6XRcjTE1jEM45\n5s+fz759+xg1alTIVYlEV1IV5y6Gmf0RuBQodc6dFuH9fOAlYEvjSy845+4LrMIQXX755aSmpmJm\nDBkyhJkzZyogJKEFfZqzNS2Ip4DHgKOt1PqGc+6KYEpKDFu3bg27BJFjFvfTnM65pcC+Fj4WzKkG\nEWkz5xylVaUc3+34wPYZxCClA842s9Vm9rKZnRLAPkXkGO07tI/M9MxA7ofRJIhByg+AQc65ajOb\nAswDTor0wenTp3/2c35+Pvn5+QEcXkTg81OcBQUFFBQUBLLPVi17b2ZDgAWRBikjfHYrcKZzbu8R\nr2vZe5EYemPbG/x8yc9ZetPSL7we6rL3ZpZn5m+Pa2bj8KGzt4VfE5GABX2KE1p3mvM5YCLQx8x2\nAL8A0gGcc78FrgJ+YGb1QDXwzUArFJFWCXqaNbQiIJxz17bw/hPAE+0porEBIiLtUFpVGugsSkiA\nmZTqs4sEo6SqhDH9xgS6T12LIdJBxGIMQgEh0kEEfSUnKCBEOoySymCv5AQFhEiHoS6GiERUWVuJ\nc46sLlmB7lcBIdIB7K7YTb/sfoFPGVBAiHQAxRXF9M/uH/h+FRAiHYACQkSiKq4opn+WAkJEIlAL\nQkSiKq5UQIhIFGpBiEhUCggRicg5p4AQkcgO1BwgxVLIzsgOfN8KCJEkF6vWAyggRJKeAkJEolJA\niEhUsZpFCQoIkaRXXFHMgO4DYrJvBYRIkovVLEpQQIgkPY1BiEhUCggRicg551eTyuoXk/0rIESS\n2N6De8lMz6RreteY7F8BIZLEYtm9AAWESFLbVbGLftmx6V6AAkIkqe0o38Hg7oNjtn8FhEgS216+\nncE9FBAiEsH2AwoIEYlCLQgRiWp7+XYG9RgUs/0rIESSVINrYOeBnQzqroAQkSOUVJbQI6NHzCZJ\ngQJCJGntOLAjpuMPoIAQSVqxHqAEBYRI0lJAiEhUCggRiUoBISJRKSBEJCoFhIhEdLDuIAdqDpDb\nLTemx1FAiCShHQd2MLD7QFIstl9hBYRIEtpRHvtJUqCAEElK8Rh/AAWESFL6uPzjmF6k1UQBIZKE\ntu7fytBeQ2N+HAWESBLasm8Lw3sNj/lxFBAiSWjz3s0M6zUs5sdRQIgkmeq6avYe3BvT+2E0aTEg\nzOyPZlZiZmuP8plfmdkmM1ttZmcEW6KINLdt/zaG9BxCakpqzI/VmhbEU8DkaG+a2VRghHPuROD7\nwG8Cqk1EIohX9wJaERDOuaXAvqN85ApgRuNnVwA9zSwvmPJE5EjxGqCEYMYgBgA7mj3fCQyM9MGG\nhgCOJtLJbdm3JW4tiLSA9mNHPHeRPpSXN51zzoHTToMLL8wnPz8/oMOLdB6b923mgqEXRH2/oKCA\ngoKCQI4VREDsAppP6RrY+NqXPPfcdB54AGbMgJwcGDsWsrICqECkE9mybwvDc6J3MfLzv/g/33vv\nvbfNxwqiizEf+DaAmY0H9jvnSiJ9cNIkeO01ePFFWL4chg6Fn/8cSiJ+WkSO1OAa/CzKnrGfRQmt\nO835HPB3YKSZ7TCzm83sFjO7BcA59zKwxcyKgN8Ct7a0z7FjYc4cePtt2LsXRo2CW26BjRvb+acR\n6eB2V+ymR0YPunXpFpfjmXMRhwuCP5CZi3asTz+FJ56AX/8azj4b7rgDzjkH7MiRDZFObunHS7nz\n9Tt56+a3Wv07ZoZzrk3fpoSYSXn88TB9OmzbBpdcAjfdBBMmwNy5UF8fdnUiiSOeZzAgQQKiSWYm\n/OAHUFgId94Jjz0GJ54IDz8MBw6EXZ1I+Dbv28ywnp00IJqkpsKVV8LSpX6s4p13YMgQ+MlPYOvW\nsKsTCU9LZzCClpAB0dy4cfDcc7BqFaSnw9e+BtOmwZtvQpyGT0QSRtHeos7bxTiawYPhwQf9OMWF\nF8L3vgdnnglPPw2HDoVdnUjsOefYsGcDJ/c5OW7HTJqAaJKVBbfeCuvXw/33w+zZcMIJfj7FrojT\ns0Q6htKqUlIshT6ZfeJ2zKQLiCYpKTBlCrz6qu9ulJf7KdzXXOPHLtT9kI4m3q0HSOKAaG7kSH/G\nY9s2OPdc+O534atfhf/5H6iqCrs6kWAUlhUysvfIuB6zQwREk+7d4fbbfffjv/8bXn7Zj138+Mf+\n1KlIMttQtkEBEYSUFLjoIpg3D1auhOxsyM+HCy6A55+H2tqwKxQ5dupixMDgwXDffbB9u7/e44kn\n/KDmv/6r75KIJIvCskJG9lELIia6dPEDmAUFsHixH5sYO9YPdP7lL1BXF3aFItHV1New88DOuK0k\n1aTTBERzo0bBI4/Ajh1w7bXw0EO+VXH33bB5c9jViXxZ0d4iTuh5Aump6XE9bqcMiCZdu8K3v+1P\ni772mp9wNX68H6uYNUsTsCRxhDH+AJ08IJo75RTfkti5E/7pn/yqVwMHwm23wfvva16FhCuMU5yg\ngPiSjAy4+mpYuNAHQ24uXHUVnH66v6r000/DrlA6ow174n+KExQQR3XCCXDPPX5c4pFH4IMP/OXn\nV17pT6HqdKnES2FZYShdjIRYUSqZHDgAf/6zv0hs/Xp/ZuTb3/ZXmWoFLIkF5xy9/rMXRT8qatN1\nGEm/olQy6d4dbr7ZX/+xYoXvglx/vT8zct99Wq9CgrerYhcZaRlxvUiriQKiHYYN812QjRvhqafg\nk0/8+hXnnOPX1ywrC7tC6QjWlKzhtNzTQjm2AiIAZn4Nzccfh+JiP59i2TIYMQKmToVnnoGKirCr\nlGS1tmQto/NGh3JsBUTA0tPh0kv9PIqdO+GGG/ziuwMH+rMhzz8P1dVhVynJZG3pWrUgOqKsLLju\nOliwwI9NTJ7sL0Hv39/P4PzLXzQZS1q2pmRNaC0IncUIQWmpv7vY3Ln+atMpU3zrYvJkv7K3SJO6\nw3V0f6A7e//fXrqmd23TPnQWI8nk5vrZmosXw4YNMHEi/OY30K+fD4rnntMy/+Jt2LOBwT0Gtzkc\n2ksBEbLcXH8Z+qJFsGWLH9R89lk/ZjF1Kvzud7p3aWcWZvcCFBAJpXdvP8fir3/1A5w33givv+6X\n1DvnHL+qd2GhrgvpTNaWhDdACQqIhNW9u5+lOXu2b0Hcc49f4GbSJB8Y//Ivfm0LrWPRsa0pDW8O\nBGiQMuk45wc2Fyzw25YtcPHF/tTq5Mn+PqfScQx+eDBLblzSrrtptWeQUgGR5IqL/eK8L7/suyMn\nn+zPikyZ4lfMSk0Nu0Jpq/2H9jPwoYEcuOsAKdb2xr7OYnRi/fv7Zf5ffNFfiv4f/wGVlX4sIy/P\nz7d4+mnYvTvsSuVYrSlZw6m5p7YrHNpLLYgObMcOv67FwoW+dTFggF/t+6KL4B/+Abp1C7tCOZqH\n3n6ILfu28PjUx9u1H3UxpEX19X4BnEWL/PbBBzBmjL/P6YUX+svVu3QJu0pp7roXruPi4Rfzna9+\np137UUDIMauq8mtxvv6634qK4Oyz4fzz/T1Exozx15VIeE567CRevOZFTs09tV37UUBIu+3d60+b\nvvEGLFniT6mefbaf5Xneeb6FkZERdpWdR9MA5f4795OWktaufSkgJHB79vgWxptv+sf1632r4txz\n/TZhAvTqFXaVHdfirYu5Z8k9LLt5Wbv3pYCQmKuogLff9mHx1lvw7rt+zc4JE3xLY/x4P4ErRefF\nAvHgWw+yu2I3D09+uN37UkBI3NXVwZo1PjT+/ndYvhz27/ddkbPO8tu4cZq41VZXP381Xx/5da4f\nfX2796WAkIRQUuLX6Vy+3Lcw3n3Xd0PGjv18GzNGXZPWGPboMF65/pVA7sWpgJCE1NDgz440hcX7\n78OqVf4K1jFj4Iwz/PbVr0LfvloVvElZdRnDfzWcfT/bF8gkKQWEJI3Dh2HTJj8PY+VKv61a5aeE\nn346jB7tt9NO8yuFH3dc2BXH38KihTzw1gMsuXFJIPtTQEhSc85fU7JqlR/XWLMG1q71NywaMgRO\nPRW+8hW/nXKKv3lRR57U9e9v/DsVtRU8eNGDgexPASEdUk2Nv6XAhx/CRx/5bd06+Phjfwbl5JN9\nK2PkSL+ddJJfUyPZuyqT/zSZW792K1eMvCKQ/SkgpFOpqfFjG4WFfn7Gxo1+6b6NG/37J57otxEj\nYPhwvw0blhzjHIcbDpPzYA6bf7Q5sBvltCcg2jdFSyQEGRmfdzmac87frKioyG+bNvkL1TZv9utm\nVFb6LsvQof7xhBP8Nniw3/r2Df/y+DUlaxiQPSCUu2hFooCQDsPMz7s4/ng/getIFRV+CvnWrb6b\nsm0bvPOOv+p1+3Y/3bxvX78e6IABn2/9+/sFhfv39+937x67lsiy7cs4d/C5sdl5GyggpNPIzvZn\nR06LsoJbbS3s2uXXA9250w+c7toF773n19PYvdvfXrG21q+1kZf3eSA1bb17f7716gU5OdCzJ3Tt\n2rpQWbZjGZeeeGmwf/B20BiEyDGqrvaTwkpK/CI9n37quzZlZf7nPXv8tm+fb5Xs3++7Pz16fL5l\nZ/uWSFaW37p1g66ZjsfSB/KjrKUM7DaM447z3amMDH/WJj3dP6al+Z9TU/3U9qbHlBQfQk0b+BZP\nTo7GIETiJjPTj2MMHdr63zl40N/rpLzcbxUVfqus9FtVFeys3IZrcFTsGMrKg34w9uBBP629psY/\n1tb6tT3q6vxEtMOH/dbQ4DfnPl/13DmYPr19f1a1IEQSxDOrn2HBxgXM/b9zA91vTNekNLPJZlZo\nZpvM7GcR3s83s3IzW9m4/bwthYh0dok2QAktdDHMLBV4HJgE7ALeNbP5zrn1R3z0DedcMLM6RDqp\n17e+zg/H/TDsMr6gpRbEOKDIObfNOVcHzAa+HuFzCT79RCSxbd23lcraylBvkhNJSwExANjR7PnO\nxteac8DZZrbazF42s1OCLFCkM1i0ZRGThk3CEmyqZ0tnMVozqvgBMMg5V21mU4B5wEntrkykE1m0\nZRGXnXhZ2GV8SUsBsQsY1Oz5IHwr4jPOuYpmP79iZr82sxzn3N4jdza92TmX/Px88vPz21CySMdy\nuOEwi7cu5pFLHglkfwUFBRQUFASyr6Oe5jSzNGADcCFQDLwDXNt8kNLM8oBS55wzs3HAXOfckAj7\n0mlOkQje3fUu33npO3x060cx2X/MLtZyztWb2W3AQiAV+INzbr2Z3dL4/m+Bq4AfmFk9UA18sy2F\niHRWi7Ys4qJhF4VdRkSaKCUSsvNnnM8dE+7g0pNicw2Gbt4rkqTKD5XzfvH7TBwyMexSIlJAiITo\nlaJXOO+E88jqkhV2KREpIERCNK9wHleOvDLsMqJSQIiEpKa+hoWbF3L5yMvDLiUqBYRISAq2FTCq\nzyj6ZvUNu5SoFBAiIXlpw0t8fWSkS5sShwJCJAQNroGXNrzElScn7vgDKCBEQrFi5wq6Z3QP5N6b\nsaSAEAnBs2uf5bpTrwu7jBZpTUqROKs7XMfcj+ay/LvLwy6lRWpBiMTZ3zb/jRE5IxjWa1jYpbRI\nASESZ39a+yduGH1D2GW0igJCJI4qaip4edPLXP2Vq8MupVUUECJx9ML6Fzhv8HkJc+/NliggROLo\nyfee5Htjvhd2Ga2mgBCJk5W7V1JcURyzdR9iQQEhEie/ee83fP/M75OWkjyzC5KnUpEkVn6onLkf\nzaXwtsKwSzkmakGIxMHM1TO5ZMQlCX3lZiRqQYjEWH1DPY+seIQZV84Iu5RjphaESIzN/Wgu/bP7\nJ9yNeVtDASESQ845Hlj2AHede1fYpbSJAkIkhv530/+SYilMGTEl7FLaRAEhEiPOOe5fej93nXtX\nwt2Ut7UUECIx8tKGl6iqreKqU64Ku5Q201kMkRiob6jnztfu5OFLHiY1JTXsctpMLQiRGPjDB39g\nQPcBTB4xOexS2kUtCJGAVdRUcO8b97Lg2gVJO/bQRC0IkYD925J/45IRl3Bm/zPDLqXd1IIQCdB7\nxe8x+8PZfHTrR2GXEgi1IEQCUt9Qz/cXfJ8HL3qQ3pm9wy4nEAoIkYA8+NaD5HTN4VujvxV2KYFR\nF0MkAMt3LufRFY/y3vfeS/qByebUghBpp/JD5Vz3wnU8eemTDOoxKOxyAmXOufgcyMzF61gi8dLg\nGrj6+avpk9mHJy97MuxyIjIznHNtataoiyHSDvcW3EtxRTHP/p9nwy4lJhQQIm0058M5zFg9gxXf\nXUFGWkbY5cSEAkKkDRZtXsTtr9zOom8tIi8rL+xyYkaDlCLHaOnHS7n+xet58ZoXOb3v6WGXE1MK\nCJFjsGTrEqbNncasabOScgm5Y6WAEGmlF9a9wDV/voY5V81h0rBJYZcTFxqDEGmBc45fvv1LHnr7\nIRbesJAz+p0Rdklxo4AQOYqq2ir+cf4/UrS3iOXfXc7gHoPDLimu1MUQiWLFzhWc8dszyEzPZOlN\nSztdOIBaECJfUl1Xzf1v3s/vV/6eJ6Y+kdRrSraXAkKkkXOO+Rvm888L/5lxA8ax6pZV9MvuF3ZZ\noVJAiABvbHuDuxffTfmhcn5/+e+5cNiFYZeUEBQQ0mk1uAYWbFjAf/39v9hduZvpE6dz3WnXJfUq\n1EFTQEins7tiN0+veprfffA7crrm8NOzf8q0U6aRlqKvw5H0NyKdQkllCQs2LmD2h7N5f/f7TBs1\njTlXzWFs/7EdaoGXoGk9COmQag/X8s6ud3hty2ss3LyQ9Z+u5+LhF3PNV65h6olT6ZreNewS46Y9\n60G0GBBmNhl4BEgFfu+c+88In/kVMAWoBr7jnFsZ4TMKCIkJ5xzb9m/jg90f8G7xu7y9823eL36f\nkX1GMmnoJCYNm8TEIRPpktol7FJDEbOAMLNUYAMwCdgFvAtc65xb3+wzU4HbnHNTzews4FHn3PgI\n+0qqgCgoKCA/Pz/sMo5JR6+5uq6arfu2UrS3iKK9RRSWFbKubB0fln5IVpcszux3JmP7j+WsAWcx\nYdAEumd0D7XeRBHLFaXGAUXOuW2NB5oNfB1Y3+wzVwAzAJxzK8ysp5nlOedK2lJQokjGfwjJWPPi\nJYs5ddyp7D24l7LqMkqrSimtKuWTyk/4pPITdlXsYteBXXxc/jGVtZUM7jGYETkjGNFrBGf2P5Mb\nRt/Aqbmnxm2Z+WT8O26PlgJiALCj2fOdwFmt+MxAIKkDorNxznHYHeZww2HqG+qpa6jzj4frqGuo\no/Zw7Wc/19TXUHO45guPh+oPcaj+EAfrD3Kw7iDVddVU11VTWVtJVV0VlbWVVNRWUFFTQXlNOeWH\nytl/aD8Vb1bweObj9M7sTZ/MPuR2y+X4zOPpm9WXUX1GcdGwixjQfQCDewwmt1suKaarA+KppYBo\nbZ/gyOZLxN+7/LnLW7m7ZjtqQ7fEtaLsSPtt/nub1mxi+Z+WR/1M0+9Hex7tMw73hWM3PY/22PSZ\nBtfw2euRfm5wDXz6zqc8+9izNLiGL22HGw77x8YQaP5Y31BPg2sgxVJIS0n7bEtPSSctJY0uqV1I\nT033jynpZKRl0CW1C11Su3Bc2nFkpGbQNb2rf0zrStf0rmSmZ5KdkU3frL5069KNrC5ZZHfJJjsj\nmx4ZPehxXA96HteThw89zL0/u7fF/14SjpbGIMYD051zkxuf3wU0NB+oNLMngQLn3OzG54XAxCO7\nGGaWPAMQIh1MrMYg3gNONLMhQDFwDXDtEZ+ZD9wGzG4MlP2Rxh/aWqCIhOeoAeGcqzez24CF+NOc\nf3DOrTezWxrf/61z7mUzm2pmRUAVcFPMqxaRuIjbRCkRST4xHxI2s8lmVmhmm8zsZ7E+XluY2SAz\nW2JmH5nZh2b2o8bXc8xskZltNLO/mVnPsGttzsxSzWylmS1ofJ7o9fY0sz+b2XozW2dmZyVBzT9p\n/Dex1sxmmVlGItVsZn80sxIzW9vstaj1mdldjd/FQjO7uKX9xzQgGidaPQ5MBk4BrjWzUbE8ZhvV\nAT9xzn0FGA/8sLHOO4FFzrmTgNcbnyeSHwPr+PysUaLX+yjwsnNuFDAaKCSBazazAcDtwJnOudPw\n3exvklg1P4X/fjUXsT4zOwU/jnhK4+/82qyF88bOuZhtwATg1WbP7wTujOUxA6p7Hn72aCGQ1/ha\nX6Aw7Nqa1TgQeA04H1jQ+Foi19sD2BLh9USueQCwHeiFH69bAFyUaDUDQ4C1Lf2dAncBP2v2uVeB\n8Ufbd6y7GJEmUQ2I8THbpfGMzRnACvxfctMZmRIgkW6h9DDwU6Ch2WuJXO9Q4FMze8rMPjCz35lZ\nNxK4ZufcLuCX+JAoxp+hW0QC19woWn398d/BJi1+H2MdEEk1AmpmWcALwI+dcxXN33M+chPiz2Nm\nlwGlzl8UF/H0cSLV2ygNGAP82jk3Bn/G6wtN80Sr2cx64S8lGIL/cmWZ2Q3NP5NoNR+pFfUdtfZY\nB8QuYFCz54P4YoIlDDNLx4fDM865eY0vl5hZ38b3+wGlYdV3hLOBK8xsK/AccIGZPUPi1gv+v/tO\n59y7jc//jA+MTxK45knAVufcHudcPfAivtucyDVD9H8HR34fBza+FlWsA+KziVZm1gU/QDI/xsc8\nZuZXDPkDsM4590izt+YDNzb+fCN+bCJ0zrm7nXODnHND8YNmi51z3yJB6wVwzn0C7DCzkxpfmgR8\nhO/XJ2TNwMfAeDPr2vhvZBJ+UDiRa4bo/w7mA980sy5mNhQ4EXjnqHuKwwDKFPwl40XAXWEO5hyl\nxnPxfflVwMrGbTKQgx8I3Aj8DegZdq0Rap8IzG/8OaHrBU7HLxmwGv9/4x5JUPN0/NXLa/FXLacn\nUs34FmQxUIsf77vpaPUBdzd+FwuBS1ravyZKiUhUunZWRKJSQIhIVAoIEYlKASEiUSkgRCQqBYSI\nRKWAEJGoFBAiEtX/B3YYHgcI0AhwAAAAAElFTkSuQmCC\n"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1116de910>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "viewer.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1106e4550>"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAESCAYAAAAFYll6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X18zfX7wPHXZTY3s7FRNHdTUXwpIje52UGySil3xTdU\nKrqRbijdfKNf36JUSimSKBWS6iuKSg2RRIiQm5rcRYxkmLHr98fnWMPGDueczznb9Xw8zuOcz+d8\nbq5zzK6970VVMcYYY3xRxO0AjDHGhB9LHsYYY3xmycMYY4zPLHkYY4zxmSUPY4wxPrPkYYwxxmeW\nPIwxxvjMkocxxhifFXU7AF+ISHvgaiAWGKuqX7ockjHGFEoSjiPMRaQM8Lyq3uZ2LMYYUxi5Xm0l\nIm+JyHYRWXHc/mQRWSMi60Tk4eNOexx4NXhRGmOMycn15AGMA5Jz7hCRCJzkkAzUArqKSE1xPAt8\nrqrLgh+qMcYYCIE2D1WdJyKJx+1uCKxX1VQAEZkEtAcuB1oDsSJyvqqODmKoxhhjvFxPHnmoCGzK\nsb0ZaKSqfYFX3AnJGGPMUaGaPE67FV9Ewq8HgDHGhABVlfweGwptHrnZAlTOsV0Zp/SRL6pqDz89\nBg0a5HoMBeVh36V9n6H88FWoJo/FQHURSRSRKOAGYJrLMRljjPFyPXmIyERgAVBDRDaJyC2qehi4\nB5gFrAImq+pqN+M0xhjzD9fbPFS1ax77Pwc+D3I45jgej8ftEAoM+y79y75Pd4XlCPOTERE9nc+0\n6a9NzFg3g3m/z+PWurfS+tzWAYjOGGNCk4igPjSYu17ycMuRrCMs3LyQGetmMH3tdLb+vZXk85O5\nNOFSen7Sk061OjGk9RBKRJZwO1RjQp5Ivn/nmBDgj0JDoSp5pB1IY9b6WcxYN4OZ62dSMbYi7aq3\n4+oaV9OoYiMiikRkH3f3Z3ez7I9lTLh+Ag0SGgTzIxgTdrx/tbodhsmHvP6tfC15FOjkoar8/OfP\nzFg7gxnrZrDsj2UkJSbRrno7rqp+FZVLVz7ptSatnES/mf24q8FdPNr8USIjIoPxEYwJO5Y8wocl\njzyIiM5YOyM7YQBcXf1q2tVohyfR43M11Na/t9JrWi927t/JhOsncGG5CwMRtjFhzZJH+PBX8nC9\nq24gDP12KFXLVGVGtxn81u83Rl49kiurX3la7RcJMQl81u0zetXrRfNxzRnx/QiyNCsAURtjAuGX\nX36hbt26xMbG8sorr3DnnXfy3//+1+/3SUxMZPbs2X6/7lVXXcWECRP8ft0zVSBLHnv2KKVL+//a\n69PW0+PjHpSILMG49uOoUrqK/29iTBgK5ZJHr169KFOmDC+88EJA71OtWjXGjh1Lq1atAnqfk/F4\nPHTv3p1evXrleYyVPE7iqacCc93z489n7i1zubza5TR4owETlk8I2f8wxhjHxo0bqVWrltthBEVQ\ne725PZ9KAOZn0XLlVFev1oBaum2p1n6ttnaY3EF37NsR2JsZE+KcXyWhp2XLlhoREaHFixfXmJgY\nXbt2rfbs2VMff/xxVVUdOnSoNmrUSA8fPqyqqq+99pr+61//0oyMDD1y5IgOGTJEzzvvPC1btqx2\n6dJF09LSsq/9zjvvaJUqVbRs2bL69NNPa2Jios6ePTvXOHr27Km9e/fWNm3aaExMjCYlJenGjRuz\n358/f742aNBAS5curZdeeqkuWLAg+72kpCR98803VVV13Lhx2rRpU+3fv7/GxcVptWrV9PPPP1dV\n1UcffTT7s5YqVUr79u2bayx5/Vt59+f/d60vB4fDA9Dhw1WvuEI1KyvX78hvDmQe0AFfDNBznj9H\np62ZFtibGRPCQjV5qKp6PB4dO3Zs9vbNN9+s//nPf1RVNSsrS1u0aKGDBw/WtWvXalxcnC5btkxV\nVV966SVt0qSJbtmyRQ8dOqS9e/fWrl27qqrqzz//rKVKldJ58+ZpRkaGPvDAA1q0aNGTJo+YmJjs\n4/v166fNmjVTVdVdu3ZpmTJl9N1339UjR47oxIkTNS4uLjtR5Yx/3LhxGhkZqW+++aZmZWXp66+/\nrgkJCXl+1txY8jhJ8jh0SLVmTdX//e+k36HfzE2dq9Veqqa9/tdL9x7cG5ybGhNCQj15HP3LXdVJ\nHkdLHqqqqampGh8frzVr1tShQ4dm769Zs+YxyWDr1q0aGRmphw8f1ieffDI7kaiqpqena1RU1EmT\nR87j9+3bpxEREbpp0yZ95513tFGjRscc36RJEx0/fnx2/DmTx/nnn3/MfUVEt2/fnutnzY2/kkeB\nbPOIjISXX4b774eDBwN/v+ZVm7O8z3IE4eJRFzN349zA39SYMCLin8fp3z/vk6tWrYrH42Hjxo3c\nfffd2ftTU1O5/vrriYuLIy4ujlq1alG0aFG2b9/Otm3bqFSpUvaxJUuWpGzZsie9f87jo6OjiY+P\nZ+vWrWzbto0qVY7tfFO1alW2bt2a67UqVKhwzH0B9u3bl6/P6k8FMnkAtGkDF10Ew4cH534xxWIY\nc+0YRlw5gm5Tu9F1ald+2/1bcG5uTIhzajnO/BEIM2bMYOHChbRu3Zr+/ftn769SpQozZ85k9+7d\n2Y/9+/eTkJDAOeecw6ZN/yx2un//fnbt2nWSz6/HHL9v3z7S0tKoWLEiCQkJbNy48ZjjN27cSMWK\nFX3+LMFsMC+wyQPghRecx+Z8LyN15trVaMcv9/xCzXI1aTCmAQO+GMCeg3uCF4Ax5gSaI/PkfL1z\n505uv/12xo4dy/jx4/n000/5/HNnMu8+ffrw6KOP8vvvvwPw559/Mm2as6xQp06dmD59OvPnz+fQ\noUM88cQTZGWdfPzXZ599ln38f/7zH5o0aULFihW58sorWbt2LRMnTuTw4cNMnjyZNWvW0K5dO58/\nZ/ny5dmwYYPP552OAp08zj0X7rwTHn44uPeNjormiaQnWHnnSvYc3MMFr17AiO9HcOjIoeAGYowB\njv2LXESyt3v37s11111HcnIy8fHxjB07lttuu43du3fTr18/rr32Wq644gpiY2Np0qQJixYtAqBW\nrVqMHDmSbt26kZCQQHx8PJUr5z3dkYjQrVs3nnzyScqWLcvSpUt59913AShbtizTp0/nhRdeoFy5\ncjz//PNMnz6d+Pj4XK9zfOki53a/fv348MMPiY+P57777jv9LywfCuQgwZyfKT0dLrwQJk6EZs3c\niWnF9hU89NVDrE9bz3OXP8d1F15ns5CaAiWUBwmGgltuuYVKlSrxVKAGofnABgnmU3Q0DBsGffvC\nkSPuxFCnfB0+//fnjLxqJINSBpE0PolFWxa5E4wxJugKYmINm+QhItEi8raIvCEi3Xw594YbICYG\nxo4NVHT5c8V5V7C091Jurnsz10++nm5Tu5G6J9XdoIwxAZdbdVO4C5tqKxHpDqSp6gwRmaSqN+Zx\nnOb2mZYtg7ZtYc0aiIsLdLSntu/QPl5Y8AIjFo2gV71ePNr8UcoUL+N2WMacFqu2Ch8FotpKRN4S\nke0isuK4/ckiskZE1onI0ebuisDRvm4+V0DVrQsdO8KgQWcYtJ+UiirFIM8gVty5grQDaVzw6gW8\n8v0rZB7JdDs0Y4w5JVdLHiLSHNgHvKOqdbz7IoBfgMuBLcAPQFegPrDbW/KYqKpd87hmriUPgF27\noGZNmD0b6tTx/+c5Ez9t/4kBXw7gt92/8Vyb52h/QfsCV8w1BZeVPMJHgVkMSkQSgU9zJI8mwCBV\nTfZuD/QeOgJ4FTgIzFPViXlcL8/kATByJEyd6iSQUPzdPGv9LPp/2Z/aZ9fm/Q7vWwIxYcGSR/jw\nV/Io6teo/CNn9RTAZqCRqu4Hbs3PBTweD4mJiSQmJuLxePB4PNnv9e4No0c7CaRTJ3+G7R9tz2+L\nJ9FDi/EtGLZgGA81fcjtkIwxBVBKSgopKSmkpqaSmprq8/mhWPLoCCSr6u3e7ZtwkkfffF7vpCUP\ngJQUuPlmWLUKvFPDhJzf//qdhmMaMqnTJDyJHrfDMeakrOQRPgpEg3ketgA5h2pWxil9+I3HA40a\nOeM/QlWV0lWYcP0Euk3txpa9W9wOxxhjjhGKyWMxUF1EEkUkCrgBmObvmwwbBiNGwHHzkYWUNue1\n4e5L76bLh11sahNjQlCRIkX49ddf3Q7DFW531Z0ILABqiMgmEblFVQ8D9wCzgFXAZFVd7e97V6kC\n/fpBjkk0Q9IjzR8hvkQ8D31pbR/GhKLCWl3navJQ1a6qmqCqxVS1sqqO8+7/XFUvUNXzVXVIoO4/\nYAAsXgxffx2oO5y5IlKEd657h0/XfsqklZPcDseYsJKYmMgLL7zAxRdfTJkyZbjxxhvJyMgAYPr0\n6dStW5e4uDiaNm3KihXOcLNx48Zx7bXXZl+jevXqdOnSJXu7cuXKLF++nBYtWgBw8cUXExMTw5Qp\nUwAYM2YM1atXp2zZsrRv355t27Zln1ukSBFGjx5NjRo1iIuL45577gn4dxAwvqwcFQ4PfFzRbOpU\n1dq1VTMzfTot6JZuW6rlniunP+/42e1QjDmBr//vgiUxMVEbNWqk27Zt07S0NK1Zs6aOGjVKf/zx\nRz377LN10aJFmpWVpW+//bYmJibqoUOHdMOGDVqmTBlVVd2yZYtWrVpVK1eurKqqGzZs0Li4uOzr\ni4hu2LAhe3v27Nlarlw5Xbp0qWZkZGjfvn21RYsWxxx/zTXX6F9//aW///67nnXWWTpz5swgfRuO\nvP6tsJUEfXP99XD22TBqlNuRnFzdCnUZ1mYYHSZ3YG/GXrfDMSZs3HvvvVSoUIG4uDiuueYali1b\nxpgxY+jduzeXXnopIkKPHj0oVqwYCxcu5NxzzyUmJoalS5cyd+5c2rZtS0JCAr/88gtz5szJLnHk\n5r333qNXr17UrVuXqKgohgwZwnfffZe9JgjAwIEDiY2NpXLlyrRs2ZJly5YF42vwu1Ac5xFUIs6S\nta1awY03QrlybkeUt5vr3syCTQvoNa0XH3T6wAYQmrAhT/rnZ1UH+d6+cPyyrVu3biUtLY23336b\nV155Jfu9zMzM7KVfk5KSSElJYf369SQlJVGmTBnmzJnDd999R1JSUp732rZtGw0aNMjejo6OpmzZ\nsmzZsiV7qdnj48m5hGw4KfTJA6B2bejaFR5/PPRLICOuHEGzt5oxfOFwHmjygNvhGJMvp/NLP5Aq\nV67MY489xqOPPprr+0lJSUybNo3U1FQee+wxypQpw7vvvsvChQvp2zfvIWcJCQnHDLhLT09n165d\np7WkbKgr9NVWRw0eDJ98AkuXuh3JyRUvWpypXaby3PznmLtxrtvhGBNW1Nsz6vbbb2fUqFEsWrQI\nVSU9PZ0ZM2ZklwKSkpL45ptvOHjwIAkJCTRr1oyZM2eSlpZGvXr1sq93/LKvXbt2Zdy4cSxfvpyM\njAweffRRGjdunF3qyCuecGTJwysuDp56ylk0KtT/PauWqcr468bTdWpXtv297dQnGGOAf9bVqF+/\nPmPGjOGee+4hPj6e6tWr884772QfV716dWJiYmjevDkAsbGxnHfeeTRt2vSY6uLBgwfTs2dP4uLi\n+PDDD2ndujVPPfUUHTt2JCEhgd9++41JkyYdc//c4glHrk9P4m/5mZ4kL0eOOCPP+/aFnj39HFgA\nPJnyJF/99hVf9/iayIhIt8MxhZhNTxI+Csysuv52JskDYMkSuPpqWLkytBvPAbI0i3bvt6NmuZq8\n0PYFt8MxhZglj/BRkOe2clX9+k6vqwED3I7k1IpIEd7t8C4frfmIKT9PcTscY0whYiWPXPz9N/zr\nX/D229CypZ8CC6AlW5eQ/F4y826Zx4XlLnQ7HFMIWckjfFjJI4BiYuCVV6BPHzh40O1oTq1+Qn2G\ntB5Ch8kd2HcoPPuMG2PCi5U8TqJDB7joIqcbbzjo9b9epGemM7HjxLDtwWHCk5U8woc1mOfBn8lj\n82aoWxe+/RYuDIPaoAOZB2j6VlN6XtyTfo37uR2OKUQseYQPSx558GfyAGfNj6lTndUHw+GP+V93\n/0qTsU34qMtHNK3S1O1wTCFhySN8WPLIg7+Tx5Ej0Lgx3HUX3HKL3y4bUDPWzqD39N4suWMJ5UuV\ndzscUwhYNWl4seSRC38nD3CmLElOdsZ+nHWWXy8dMI/Nfoxl25cxvet0+49tjDkl620VAPXqQffu\n8OCDbkeSf4M9g/lj3x+M+XGM26EYYwqgsCp5iEh74GogFhirql/mcozfSx4A+/Y5Yz/GjoXLL/f7\n5QNi1Z+raDGuBd/f9j3nxZ/ndjjGmBBWKKqtRKQM8Lyq3pbLewFJHgDTp8P998NPP0GJEgG5hd8N\n/244H67+kLk3zyWiSITb4RhjQlRYVFuJyFsisl1EVhy3P1lE1ojIOhF5+CSXeBx4NbBRnqhdO6fr\n7jPPBPvOp69f434UiyjGsAXD3A7FGFOAuFLyEJHmwD7gHVWt490XAfwCXA5sAX4AugINgEuAYcA2\nYCjwharOzuPaASt5AGzdChdfDHPmQK1aAbuNX/3+1+/Uf6M+X3b/kroV6rodjjEmBIVFyUNV5wG7\nj9vdEFivqqmqmglMAtqr6gRVvV9VtwJ9gdZAJxHpHdyoHQkJ8OST0Ls3ZGW5EYHvqpSuwvNtnqf7\nx93JOJzhdjjGmALAtTYPEUkEPs1R8ugEtFXV273bNwGNVDXvNR9zv64mJSWRmJhIYmIiHo8Hj8fj\n19iPHIGmTeG225xHOFBVOn7QkfPjz+e5Ns+5HY4xxmUpKSmkpKSQmppKamoqc+bMCY8G81ySR0cg\n2R/JIxif6aefnF5XK1ZA+TAZh/dn+p9cPOpiJnWaRIuqLdwOxxgTQsKi2ioPW4DKObYrA5tdiuWU\nLrrIGXH+wANuR5J/Z0Wfxeh2o7n5k5v5O+Nvt8MxxoSxUEoei4HqIpIoIlHADcA0l2M6qSeegO++\ng1mz3I4k/6654BpaVWvF/bPudzsUY0wYc6ur7kRgAVBDRDaJyC2qehi4B5gFrAImq+pqN+LLr+ho\nGDnSmfdq/363o8m/4W2H8/VvXzPtl5DOzcaYEBaWgwRPJlhtHjndeCOce254jf+Yt3EeN3x4A8v7\nLOes6DCZsMsYEzCFYoT5ybiRPP74w2kDmT0b6tQJ6q3PyENfPsT6tPVM7TLVJk80ppAL5wbzsFWh\nAjz1VHiN/QB4quVTrE9bzzvL33E7FGNMmLHk4Se33+4sFjUmjCaxLVa0GBOun0D/L/uzcc9Gt8Mx\nxoQRq7byo5UroVUrWL4czjnHlRBOy9Bvh/LFhi/4qsdXFBH7e8KYwsiqrVxUu7Yz4vzuu+HwYbej\nyb8Blw3g0JFDvLzwZbdDMcaECSt5+NmBA3Dddc7zxIlQsaJrofhkQ9oGGo9tzJyb51DrrDCZ8dEY\n4zdW8nBZiRLw+efQti00aBA+AwjPiz+PZ1o9Q/ePu3PoyCG3wzHGhDgreQTQnDnw739Dz57OTLxF\ni7od0cmpKtdMvIZ6FerxVKun3A7HGBNENs4jhJIHwI4dcNNNkJHhVGMlJLgd0cn9se8P6o6qyyc3\nfkLjSo3dDscYEyRWbRVizj7bqcZq0wbq14cvvnA7opOrUKoCI68aSY+Pe5B+KN3tcIwxIcpKHkGU\nkuJUY91yCwweHNrVWN0/7k5sVCwjrx7pdijGmCCwaqsQTh4A27c71ViZmfD++6FbjbXn4B5qv1ab\nj2/4mEsrXup2OMaYALNqqxBXvjzMnOkMJqxfH7780u2IclemeBkebvowQ74d4nYoxpgQZCUPF339\nNXTvDr16waBBEBHhdkTH2p+5n3NfPpfZPWbzr7P/5XY4xpgAspJHGGnVCpYsgQULnCVtt21zO6Jj\nlYwsyX2N72Po/KFuh2KMCTGWPFxWoYIzkLBly9CsxrqzwZ18vu5zft39q9uhGGNCiFVbhZCvv3Ya\n02+7LbSqsf7z9X/4c/+fjGo3yu1QjDEBUuCrrUQkWkR+EJGr3Y7F31q1gh9/hG+/daqxtm51OyJH\nv8b9+ODnD9iyd4vboRhjQkTYJQ/gIWCy20EESoUKTtVVq1ZwySXOAEO3lStZjpvr3syL373odijG\nmBDhSvIQkbdEZLuIrDhuf7KIrBGRdSLycC7ntQFWAX8GK1Y3RETAf/4DH3wAd9wB/fvDIZfnKnyw\nyYOMXz6enft3uhuIMSYkuFXyGAck59whIhHAq979tYCuIlJTRLqLyHARSQCSgMZAN+B2KeALb7do\nAcuWwS+/QLNm8KuLbdYVYyvSqWYnRnw/wr0gjDEhw7UGcxFJBD5V1Tre7SbAIFVN9m4PBFDVE/qJ\nikhP4E9V/SyX98K2wTwvqjBiBDz9NLz6KnTp4k4cR9f82HDvBmKLxboThDEmIHxtMA+l2ZUqApty\nbG8GGuV2oKq+fbILeTweEhMTSUxMxOPx4PF4/BelC0SgXz+n9HHjjfDVV/DSS1CyZHDjOC/+PNqe\n15bXf3idh5udUKtojAkjKSkppKSkkJqaSmpqqs/nh1LJoyOQrKq3e7dvAhqpal8fr1vgSh457d0L\nd97prJM+eTL8K8gDv1fuWEmbCW349d5fKRFZIrg3N8YETDh31d0CVM6xXRmn9GFyiI2Fd9+FBx4A\njwfefNOp1gqW2mfXpnGlxoxdOjZ4NzXGhJxQSh6LgeoikigiUcANwDSXYwpJInDrrTB3rtMW0rWr\nUyIJlkeaPcKwBcNsuVpjCjG3uupOBBYANURkk4jcoqqHgXuAWTjdcSer6mo34gsXNWvC999DXJwz\nJmTx4uDct2HFhtQoW4P3fnovODc0xoQcm56kgJgyBe6+GwYOhPvvd0ongZSSmsIdn97B6rtXE1Ek\nROZRMcactnBu8zBnoHNnpxTywQdwzTWwM8Bj+ZKqJlGuZDmmrp4a2BsZY0KSJY8CpFo1mDfP6YFV\nr57TJhIoIsJjzR/jmXnPUBhLesYUdpY8CpjISHj2WRgzxhlMOGxY4HpjXVX9KgA+W3fCWE1jTAFn\nyaOASk6GRYvgww+hQwfYs8f/9xARHm3+KE/Pe9pKH8YUMnkmDxHZJyJ/5/EIYsdQc7qqVHGqsSpV\nggYNnHmy/K1jzY7sOrCLORvn+P/ixpiQZb2tColJk6BvXxg61Fkz3Z/GLR3H+yvf58vuIbYMojEm\n33ztbZXv5CEiZwPFj26r6u++hxd4ljzytno1dOwIjRvDyJFQwk+zixw6cojqr1RnSucpNKzY0D8X\nNcYEld+76orItSKyDvgNmAOkAiGwRJHxVc2aTjvIwYPQpAmsX++f60ZFRDHgsgEM+XaIfy5ojAl5\n+Wkw/y/QBFirqtWA1sD3AY3KBEypUvDee84iU5ddBh9/7J/r9qrXi4WbF7Jyx0r/XNAYE9Lykzwy\nVXUnUEREIlT1G6BBgOMyASQCd90F06c7o9EHDIDMzDO7ZonIEtzX6D6GfnvC8ivGmAIoP8ljt4jE\nAPOA90RkBLAvsGGZYGjYEJYsgZUroXVr2Lr1zK5356V3MmvDLDakbfBPgMaYkJWf5HEdsB+4H5gJ\nrAeuCWRQJnjKloUZM6BNG6c7b0rK6V8rtlgsdza4k+fmP+e3+Iwxocm66ppsX30F3bs7qxY+9BAU\nOY0hpDv37+SCVy/gpz4/UTG2ov+DNMYEhN+76orIPuDoQVFAJLBPVUNyEWtLHmdm82ZnWpOyZeGd\nd5zp3n314KwHUZQX277o/wCNMQHh9666qlpKVWNUNQYoAXQAXjuDGE0Iq1TJqbo67zyoXx9+/NH3\nazx42YOMXzaenfsDPLWvMcY1PlVMqGqWqn4CJAcoHhMCoqLgpZec0eht28Lbb/t2fkJMAl3+1YWX\nF74cmACNMa7LT7VVxxybRYD6QJKqNglkYKfLqq38a9UqZ2LFli2dhFKsWP7O+3X3rzQc05AN926g\ndPHSgQ3SGHPGArEY1DVAO+/jCuBvoP3phXdmxPG0iIwQkR5uxFDY1KrljErfvh2Skpw2kfw4N+5c\nrqx+Ja8vfj2wARpjXBFWva1E5HqcxLUT+ExVv87lGCt5BIAqPPecU/p47z1o1erU56z6cxWe8R4W\n37GYKqWrBD5IY8xp81tvKxF5JcemApLjNap67xkE+RZwNbBDVevk2J8MvAREAG+q6rPHnfcwkKaq\nY0Rkiqp2zuXaljwCaPZsuOkmeOAB6N//1GulD5s/jI/WfMScm+cQFREVnCCNMT7zZ7XVEu+jGHAJ\nsBZYB9TD6bJ7JsZxXKO7iEQAr3r31wK6ikhNEekuIsNFJAHYDBxd1ijrDGMwp6F1a2et9ClTnHXT\n//775Mc/eNmDlCtZjoFfDQxOgMaYoMhPg/n3QDNVzfRuRwLfqmqjM7qxSCLw6dGSh4g0AQaparJ3\neyCAqg7NcU4J4BWcEe+rVfWECnUreQRHRoYzmHDOHPjoI2fG3rykHUjjktGXMLztcK6veX3wgjTG\n5JuvJY+i+TimDBAL7PJux3j3+VtFYFOO7c3AMQlKVQ8At53qQh6Ph8TERBITE/F4PHg8Hr8Gapxe\nV6NGwbhx0KIFvP46dOqU+7HxJeKZ3Gky10y8hosrXMy5cecGN1hjzAlSUlJISUkhNTWV1NRUn8/P\nT8njFmAwkOLdlQQMVtXxPt/t2OsmcmzJoyOQrKq3e7dvAhqpal8fr2sljyBbssRJHJ07wzPPQNE8\n/iR5eeHLTPhpAvNvnU+xovns82uMCYpAjDAfBzQGPgY+AhqfaeLIwxagco7tyjilDxPi6teHxYth\n+XJngsUdO3I/7t5G91K1TFUe/OLB4AZojPG7PJOHiNT0PtcHzsGpUtoMJIjIJQGIZTFQXUQSRSQK\nuAGYFoD7mAAoWxY++wyaNXNm51248MRjRIS3rn2Lz9d/zuSVk4MfpDHGb07WVXeMqt4uIin8MzFi\nNlVtedo3FZmIU/1VFtgBPKGq40TkSv7pqjtWVX1e19Sqrdw3bRrcdhs8+ST06XNid94ft/1I23fb\nMv/W+dQoW8OdII0xx/D7rLrhxpJHaFi3zpnWpH59pzG9RIlj33/9h9cZtWQUC3stpERkidwvYowJ\nGr+3eYhkwwgEAAAcGklEQVRIZxGJ9b7+j4h8FKBqK1OAVK/uVF1lZDhVWRs3Hvt+nwZ9qHVWLfrN\n7OdOgMaYM5Kfua2eUNW9ItIMaA28BYwKbFimIIiOhvffh3//Gxo3hq9zTCYjIrzR7g3mbJzDuz+9\n616QxpjTkp/kccT73A4Yo6rTcRaEMuaURJypTI4mkRdecObJAogpFsOUzlO4f9b9rPpzlbuBGmN8\nkp/ksUVE3sDp/TRDRIrn8zxjsrVs6UxrMnEidOsG6enO/ovKX8TQ1kPpPKUz6YfS3Q3SGJNv+UkC\nXYCZwBWqugeIAwYENCpTIFWpAvPmOaPTmzSBDRuc/bfWu5UGCQ2467O7sM4OxoSH/AwSTAf+BJp5\ndx0G1gcyKFNwlSjhTGnSuzdcdhnMnOm0f7x21Wss3rqYccvGuR2iMSYf8jM9yWCc1QMvUNUaIlIR\n+EBVmwYhPp9ZV93w8e23cMMNcPfd8MgjsHrnKpLGJzG7x2wuKn+R2+EZU6gEYiXBowswpQOo6hac\nyRGNOSPNmjmrFH76KXTsCJWK1eLFK16k85TO/J1xirnejTGuyk/yyFDV7LUzRCQ6gPGYQqZiRUhJ\ngbPPhkaNoGHx7rSo0oI7pt9h7R/GhLCTJg8REWC6iIwGyojIHcBs4M1gBGcKh6PTuz/4IDRvDm2O\njGDVn6sYvWS026EZY/Jw0jYPb/JYAdwPtPXunqWqXwYhttNibR7h7fvvnandr7l5LR/ENGXWTbO4\n5Byb0MCYQPP73FYi8jYwUlUXnWlwwWDJI/xt3+5d4rbqZPbUf5RlfX6kdPHSbodlTIEWiOTxC3A+\nsBFvozmgqhqS3WEseRQMmZnQvz+M//Mu2rTN5MOeY9wOyZgCLRDJIzG3/aqa6ktgwWLJo2B5fXwa\nd/9yPiPrLOfObpVPfYIx5rTYlOyWPAqcmyY8yLRpyp3nvsgzz0BEhNsRGVPwWPKw5FHgbN67mTqv\nXcRFKesollWWiROdlQuNMf4TiEGCxriqUmwlOtS8npYPj+Tii51lbpctczsqYwq3sCp5iEgV4GUg\nDVirqs/mcoyVPAqgNTvX0GJcC37r9xuffhRN377w0kvONO/GmDNX0EsetYEPVbUXUM/tYEzwXFju\nQppXbc7YpWO58UaYPRsGDYL773d6ZhljgsuV5CEib4nIdhFZcdz+ZBFZIyLrROThXE79HuglIrNx\npok3hcjDTR/m+QXPk3kkk4sucubFWrMGrrgCduxwOzpjChe3Sh7jgOScO0QkAnjVu78W0FVEaopI\ndxEZLiIJwM3AIFVtDVwd5JiNyxpWbEj1stWZuHIiAPHxMH26M7X7pZfC4sUuB2hMIeJK8lDVecDu\n43Y3BNaraqqqZgKTgPaqOkFV71fVrTiljXtF5HXgt+BGbULBwKYDeXb+s2R55+qMiICnn4bhw+HK\nK2H8eHfjM6awKOp2ADlUBDbl2N4MNMp5gKr+DHQ+1YU8Hg+JiYkkJibi8XjweDx+DdS45/JzL6d4\n0eJMXzuday+4Nnt/hw5w4YVw/fXwww9OMomKcjFQY0JcSkoKKSkppKamkpqa6vP5rvW28o5c/1RV\n63i3OwLJqnq7d/smoJGq9vXxutbbqoCb8vMUXlz4IgtuXYAzd+c//voLuneHtDT48EOoUMGlII0J\nM+Hc22oLkHP+ico4pQ9jjtGhZgd27t/JvN/nnfBe6dLwySfQpo0zHuS771wI0JhCIJSSx2Kguogk\nikgUcAMwzeWYTAiKKBLBQ5c9xNBvh+b6fpEiTjfe11+H9u3hjTeCHKAxhYBbXXUnAguAGiKySURu\nUdXDwD3ALGAVMFlVV7sRnwl9PS7uwbI/lrH8j+V5HnPNNc466S+/DLffDhkZQQzQmAIurEaY54e1\neRQew+YPY+kfS3m/4/snPe7vv6FnT9i2DaZOhYSEIAVoTBgJ5zYPY3zSu0FvvtjwBb/u/vWkx8XE\nOI3n11zjjAeZPz9IARpTgFnyMGErtlgsvev35vkFz5/y2CJF4NFHYcwYp1vvqFFgBVRjTp9VW5mw\ntn3fdmqOrMnqu1dTvlT5fJ2zfj1cdx00bgwjR0KxYgEO0pgwYNVWplApX6o8XWt3ZcT3I/J9zvnn\nO1149+yBpCTYsiWAARpTQFnyMGGv/2X9Gb1kNHsz9ub7nJgYmDLFKYFceinMO3HIiDHmJCx5mLBX\nLa4abc9vy+jFo306TwQGDoS33oJOneC116wdxJj8sjYPUyD8tP0nkt9N5td+v1K8aHGfz9+w4Z9S\nyGuvQXHfL2FMWLM2D1MoXVT+IuqdU48Jyyec1vnnnee0g6SnQ4sWsNkmxjHmpCx5mAJjYNOBPLfg\nOY5kHTmt80uVgkmTnCqshg1h7lw/B2hMAWLJwxQYzao04+zos/lo9UenfQ0ReOghZ12Qzp3h1Vet\nHcSY3FjyMAWGiDCw6UCGzh/KmbZ7XXGFU431xhtw661w8KCfgjSmgLDkYQqUq2tcTcbhDL769asz\nvta55zoJ5MABaN4cNm069TnGFBaWPEyBUkSK8HDThxk6P/fp2n0VHQ0TJ0KXLk47SEqKXy5rTNiz\n5GEKnBtr38j6tPUs2rLIL9cTgQEDYMIEuPFGZ4p3awcxhZ2N8zAF0ivfv0LKxhSmdpnq1+umpjrr\npNeuDaNHQ8mSfr28Ma6xcR7GAL0u6cW3v3/Lmp1r/HrdxERnSndVaNbMSSbGFEYhmzxEpJqIvCki\nU7zb0SLytoi8ISLd3I7PhLaSkSW559J7GDZ/mP+vXdKpwurRw5mZ96szb5s3JuyEfLWViExR1c4i\n0h1IU9UZIjJJVW/M43irtjIApB1I4/wR5/PTnT9RKbZSQO7xzTfQrRv07w8PPOC0jxgTjkKu2kpE\n3hKR7SKy4rj9ySKyRkTWicjD+bhUReBoZ8nTG0JsCpX4EvHcUvcWhn83PGD3aNkSvv/e6ZHVrZsz\nvYkxhUEwqq3GAck5d4hIBPCqd38toKuI1BSR7iIyXERyW2V6M1DZ+zpkq9tMaLm/yf2MWzaOtANp\nAbtHlSrOlO5RUdCkCfx68lVxjSkQAv5LWFXnAbuP290QWK+qqaqaCUwC2qvqBFW9X1W3iki8iIwC\n6nlLJh8BHUXkNWBaoOM2BUOl2Ep0rd2Vu2bcdcajzk+mRAlnSpM77nASyKxZAbuVMSEhKG0eIpII\nfKqqdbzbnYC2qnq7d/smoJGq9vXDvazNwxzjQOYBksYn0bFmRx5ulp8a0jMzbx7ccAPcey88/LC1\ng5jw4GubR9FABnMSAf3t7vF4SExMJDExEY/Hg8fjCeTtTIgrEVmCj2/4mIZvNqRO+TpcVf2qgN6v\neXNYtAg6doQlS2DcOGfGXmNCSUpKCikpKaSmppJ6Gn3O3Sp5NAYGq2qyd/sRIEtVn/XDvazkYXL1\n3abvaD+pPfNumccF5S4I+P0OHoR77oGFC+Hjj6F69YDf0pjTFnK9rfKwGKguIokiEgXcgLVjmABr\nUrkJQy8fyrWTrmXPwT0Bv1/x4jBmDPTtC02bwowZAb+lMUETjK66E4EFQA0R2SQit6jqYeAeYBaw\nCpisqqsDHYsxt9a7lbbnteXfH/37tBeN8oUI9O4Nn3ziNKb/97+QlRXw2xoTcCE/SNBXVm1lTiXz\nSCZt321Lo4qNGHL5kKDdd+tWpx3knHPg7bchJiZotzbmlMKl2soY10RGRPJB5w+Y/PNkJq6YGLT7\nJiQ4U7qfdRY0agRr1wbt1sb4nSUPUyiVK1mOT278hH4z+/Hjth+Ddt9ixZzZeO+7z5lYcfr0oN3a\nGL+yaitTqE1dNZUHvniAH27/gbOjzw7qvb/7zlknvXdveOwxKGJ/yhkX+VptZcnDFHpPfPMEKakp\nfNXjK6IiooJ6761boVMnKF/eaQeJjQ3q7Y3JZm0exvhosGcwcSXi6Pd5v6Df+2g7SPnyTjvIL78E\nPQRjToslD1PoFZEiTLh+AnN/n8uoxaOCfv+oKBg1ypnSvXlz+PTToIdgjM+s2soYr/Vp62n6VlM+\n7Pwhzas2dyWGhQudaqw77oDHH7d2EBM81uZhycOcgVnrZ3Hz/27m+9u+p0rpKq7EsG2bk0DOOgve\necfaQUxwWJuHMWeg7flt6d+kP9dNuo79mftdieGcc5wVChMSoGFDWOPfZdiN8QsreRhzHFWlxyc9\nOJx1mPc7vI+4OKf62LHwyCPw5ptw7bWuhWEKASt5GHOGRIQ32r3B+rT1PDf/OVdj6dULpk2Du++G\nJ5+0ebFM6LCShzF52Lx3M43ebMSYa8YEfA2QU/njD6cdJD4eJkyA0qVdDccUQFbyMMZPKsVWYkrn\nKdz8yc38stPdARgVKsDXX0Plyk47yGqbg9q4zJKHMSdxWeXLGNJ6CNdOupb0Q+muxhIVBSNHwsCB\n0KKFs8CUMW6xaitj8qHLlC7UPrs2TyQ94XYoAPzwgzO9e48eTltIRITbEZlwZ+M8LHmYAPht9280\nGNOAlXeu5JyYc9wOB4AdO6BLFyhZEt57D+Li3I7IhDNr8zAmAKrFVePWurfyxDehUfIAOPts+PJL\nqFEDLr0UVq50OyJTmIR08hCRaiLypohM8W63F5E3RGSSiLRxOz5TuDzW4jGmrZ3Giu0r3A4lW2Qk\nvPQSDB4MLVvCBx+4HZEpLMKi2kpEpqhq5xzbZYDnVfW2XI61aisTMK98/wrT101n1k2z3A7lBEuX\nQocOTlXW009D0aJuR2TCSUhWW4nIWyKyXURWHLc/WUTWiMg6EXnYh0s+Drzq3yiNObU+DfqQuieV\nmetnuh3KCerVcxrSlyyBK6+EXbvcjsgUZMGqthoHJOfcISIROAkgGagFdBWRmiLSXUSGi0jC8RcR\nx7PA56q6LBiBG5NTZEQkz13+HP2/6M/hrMNuh3OCcuVg5kwnkTRoAMvsf4kJkKAkD1WdB+w+bndD\nYL2qpqpqJjAJaK+qE1T1flXdKiLxIjIKqCsiA4F7gNZAJxHpHYzYjTnetRdcS7mS5Xhr6Vtuh5Kr\nokXhuedg6FBo08bpiWWMv7lZK1oR2JRjezPQKOcBqpoG9DnuvFcCHJcxJyUivHDFC7Sb2I6utbsS\nUyzG7ZBydcMNULOm0w6yeDEMG2btIMZ/3PxRClirtsfjITExkcTERDweDx6PJ1C3MoVU/YT6XH7u\n5Tw7/1n+2+q/boeTp4sugkWL4N//dkohkyc7XXyNSUlJISUlhdTUVFJTU30+P2i9rUQkEfhUVet4\ntxsDg1U12bv9CJClqs+e4X2st5UJik1/baLu6Los77OcSrGV3A7npI4cgSeegHffhalTnfYQY3IK\nyd5WeVgMVBeRRBGJAm4AprkYjzE+qVy6Mn3q9+Gxrx9zO5RTiohwuu8OH+70xHrnHbcjMuEuWF11\nJwILgBoisklEblHVwzgN4LOAVcBkVbW5Qk1YGdhsIF9s+IIft/3odij50qGDs0rh//0f3HcfZGa6\nHZEJV2ExSNAXVm1lgm304tFM+nkSX/f42tVVB32xezd06wYHDzqj0s86y+2IjNvCqdrKmAKh1yW9\n2JG+g0/Xfup2KPkWFwfTp0Pjxs68WD+GR8HJhBBLHsacoaJFivJ8m+cZ8OUAMo+ETz1QRAQMGeJ0\n4W3b1saDGN9Y8jDGD5LPT6Zq6aqMXjLa7VB81rkzzJ7t9MZ68EE4HHoD500IsjYPY/zkp+0/0WZC\nG3655xfKFC/jdjg+S0uDG2+ErCyYNMmZ6sQUHtbmYYxLLip/Ee2qt+OZec+4HcppiY+Hzz6DSy5x\n2kGWL3c7IhPKrORhjB9t/XsrdV6vw5I7lpBYJtHtcE7bpEnQty+88opTGjEFny1Da8nDuOzJlCdZ\ns2sNEztOdDuUM7JsGVx/vdMmMmSIrZNe0FnysORhXJZ+KJ0LXr2AqV2m0qhSo1OfEMJ27nQmWIyI\ncEoj8fFuR2QCxdo8jHFZdFQ0T7V8ige+eIBw/0OmXDmYNQvq1HHaQVaEzgq8xmWWPIwJgB4X9yD9\nUDpTV091O5QzVrQovPCCM6VJq1a2TrpxWLWVMQEy+9fZ3DH9DlbdtYpiRYu5HY5f/Pij0wZSt66T\nUBIT3Y7I+ItVWxkTIlqf25qa5Woy8oeRbofiN5dcAitXOsvc1q8PgwbB/v1uR2XcYMnDmAAa1mYY\nQ74dwq79u9wOxW9KlIDHH3d6Y61d66xW+MEHYAX+wsWqrYwJsDun30mxosV4Kfklt0MJiLlz4d57\noUwZGDHCWb3QhB+rtjImxDzZ8kne/eld1u1a53YoAdGiBSxZ4gwmbNMG7r4bdhWcgpbJg5U8jAmC\nIfOGMHX1VFpVa0Vc8TjiSsTl+lymeBkiioTeaLzMI5mkZ6azP3M/GYczKF+qPCUjS55wXFqaM8Hi\nBx/A4MFwxx1Ob62C4nDWYbb9vY3Nezezae8mdqTvoETREsQWi6V08dLEFos95lEqqhRFJDz+RrdB\ngpY8TAjKOJzB+yveZ3v6dnYf2M2eg3vYfXC38zjwz/PejL1ER0VTpniZY5OL93XJyJIIgoggOP/P\nc3t9dFGq3F5nHM7ITgTph9LZf3j/P68z9x/7nnc7S7OIjowmOiqayCKR7EjfQcnIklSKrUTl0pWp\nFFOJSrGVsrfTt1Xi+UGV+GtHKUaMAI/Hla/dJ4ezDvPHvj/Y9Nem7ORwzPNfTrIoV7Kc85ljK1E+\nujwHDx9kb8beYx5/ZfzF3oy97M/cT6moUickldhisZQu5iSb6MhoIiMiKVqk6AmPyCIn7s/t2HPj\nzqVK6Spn9PkLTPIQkWrAY0BpVe3s3RcNpACDVXVGHudZ8jBhK0uz2Jux95iEkvP5QOYBFM0efJjb\na8W7nctrVaVY0WJER0ZTMrIk0VHeZ+92XvuiIqKOWSVRVdl1YFf2L9XNezc7j783H7NPsqLI3FWJ\nskUr46lfiQvPcRJMhVIViCkWQ3RkNKWiShEd5X32/iI9U0eyjvBXxl/sObiHvw46zzkff2X8xe4D\nu9m2b1t2cti+bztlS5alcqyTGCrHVs5OEkf3JcQk+BTfkawj7Du074SkkvOx79A+DmcdPuaReSTz\nn23NY//RfVmZ9KrXix4X9zij76zAJI+jRGRKjuTxJPA3sNqShzGhTVXZfXA367ZvZsTbm/lk9mYa\ntN5Eldqb2XnwD9IPpbPv0D7SM73P3u2IIhG5JpXs7UjnWVXZk3FcUvAmivTM9Oy/7ssUL5P9KF28\nNGWK/fP6nFLnZCeIhJgEoiKi3P7aXBNyyUNE3gKuBnaoap0c+5OBl4AI4E1VfTaP86eoamcRaQPE\nA8WBnZY8giMlJQVPONQ5hIHC/l1u3AgDBsCiRfDQQ1C+PJQsCdHR/zxKllSKFj+ERu7jsKSTnrnv\nhCRz9LHhxw1c2vTSY5NDsdKUkDJoRgz7/i7C3r3k+fj7b6fbcdmyzpxdZcse+zo+Hor5eWznoUMn\nxpGR4cRRsmTuz8GakNLX5BGMpqxxwCvAO0d3iEgE8CpwObAF+EFEpgENgEuAYaq69bjrJAHRQC3g\ngIh8Zlki8Ar7Lzx/KuzfZdWqTkP6N9/A+PGwbx+kpx/72L9fSE8vRnp6MY4cKZudXI5PMtHRsGbN\nYFZ80IO9e+Gvv/75ZZyVBaVLQ2zsyR9nnQUHDsCff8KaNU5j/65dxz5HRR2bTI5PMLGxzjVOlaSO\nvj5yxIktJuafOKKi4OBB5zr795/4HBl58uRSsiTcdJMzA3IwBTx5qOo8EUk8bndDYL2qpgKIyCSg\nvaoOBSZ498UDzwB1ReRhVX3cu78n8KclDmPCU8uWzuNUDh8+mlBOTDLp6c4sv336nJgUihUDyfff\nz3lTdRLc8Qnl6POmTU7Sio527hsf70zXkjMxnGlsqk5pJbekkvN1jRpn/nl95VYnuorAphzbm4Fj\n5q5W1TSgz/EnqurbgQ3NGBMKihZ1/kovXTr395ctC2wvLhEnEcTEuDeHl4iTcIoVg7g4d2LIS1Aa\nzL0lj0+PtnmISEcgWVVv927fBDRS1b5+uJeVSIwx5jSEWptHbrYAlXNsV8YpfZwxXz68McaY0+PW\n0MfFQHURSRSRKOAGYJpLsRhjjPFRwJOHiEwEFgA1RGSTiNyiqoeBe4BZwCpgsqquDnQsxhhj/CPk\nBwnmV37HjZj8EZFUYC9wBMhU1YbuRhRechvf5O1BOBmoCqQCXVR1j2tBhpE8vs/BwG3An97DHlHV\nme5EGD5EpDLO0ImzAQXeUNURvv58hseMXaeQY9xIMs44kK4iUtPdqMKeAh5VrWeJ47SMw/l5zGkg\n8KWq1gBme7dN/uT2fSrwovdntJ4ljnzLBO5X1X8BjYG7vb8vffr5LBDJgxzjRlQ1E5gEtHc5poLA\nOh+cJlWdB+w+bve1wNGu5m8D1wU1qDCWx/cJ9jPqM1X9Q1WXeV/vA1bjDJ/w6eezoCSP3MaNVHQp\nloJCga9EZLGI3O52MAVEeVXd7n29HSjvZjAFRF8RWS4iY0WkjNvBhBvvMIp6wPf4+PNZUJJHwWi4\nCS1NVbUecCVOsba52wEVJN4ZEuzn9sy8DlQD6gLbgBfcDSe8iEgpYCrQT1X/zvlefn4+C0ryCNi4\nkcJKVbd5n/8EPsapGjRnZruIVAAQkXOAHS7HE9ZUdYd6AW9iP6P5JiKROIljgqp+4t3t089nQUke\nNm7Ej0SkpIjEeF9HA1cAK9yNqkCYBvT0vu4JfHKSY80peH/BHXU99jOaL+IszDIWWKWqL+V4y6ef\nz4LUVfdK/umqO1ZVh7gcUtjyLsT1sXezKPCefZ++8Y5vSgLK4dQfPwH8D/gAqIJ11fVJLt/nIMCD\nU2WlwG9A7xx19iYPItIMmAv8xD9VU48Ai/Dh57PAJA9jjDHBU1CqrYwxxgSRJQ9jjDE+s+RhjDHG\nZ5Y8jDHG+MyShzHGGJ9Z8jDGGOMzSx7GGGN8ZsnDGGOMzyx5GOMDEYkWkRkiskxEVojIQyIy1fte\nexHZLyJFRaS4iGzw7j9PRD73zlA8V0Qu8O4/S0Q+FJFF3sdl3v2DRWSCiCwQkbUicpt7n9iY3BV1\nOwBjwkwysEVVrwYQkVigt/e95jjzKzUEIoGF3v1v4EydsV5EGgGvAa2Bl4HhqjpfRKoAM3EWMwOo\njbNQTylgqYjMODpZpTGhwJKHMb75CXheRIYC01X1WxHZICIXApcCLwItcOZYm+edWPIyYIozHx0A\nUd7ny4GaOfbHeI9X4H+qmgFkiMg3OAnpf4H/eMbkjyUPY3ygqutEpB7Oetr/FZHZOJPMXYWzvOds\nnFXYigD9cZLIbu/aKMcToJGqHjpmp+S6OF6W3z6EMX5gbR7G+MA7DfhBVX0PeB64BJgH3AcsUNWd\nQFmghqr+rKp7gd9EpJP3fBGRi7yX+wK4N8e16x59CbQXkWIiUhZn9tgfAv/pjMk/K3kY45s6wDAR\nycIpafTBWQP6bJwSCMByjl3C89/A6yLyOE5byESc6q97gZEishzn/+Ic4C6caqufgG9wpiD/P1X9\nI8Cfyxif2JTsxoQYERkE7FNVW1bVhCyrtjImNNlfdSakWcnDGGOMz6zkYYwxxmeWPIwxxvjMkocx\nxhifWfIwxhjjM0sexhhjfGbJwxhjjM/+H+/5mhR/yHmkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11167a1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(fixedPoint[...,0], fixedPoint[..., 1], label=\"fixed point\")\n",
    "plt.semilogy(newton[...,0], newton[..., 1], label=\"newton\")\n",
    "plt.ylabel(\"residual\")\n",
    "plt.xlabel(\"sweep\")\n",
    "plt.legend()"
   ]
  },
  {
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}