Exampled of using solve_bvp

solve_bvp_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples of using `solve_bvp`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examples are taken from the papers:\n",
    "\n",
    "[1] \"A BVP Solver Based on Residual Control and the Maltab PSE\". \n",
    "\n",
    "[2] \"Solving Boundary Value Problems for Ordinary Differential Equations in Matlab with bvp4c\"\n",
    "\n",
    "I use same initial meshes and guesses for function values and parameters. In bvp4c both relative and absolute tolerances are used (in a somewhat confusing way) and also their RMS residuals are not normalized by the meash interval lenghts. Thus tolerance settings don't really match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import solve_bvp\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Measles spread model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equations:\n",
    "$$\n",
    "y_1' = \\mu - \\beta(t) y_1 y_3 \\\\\n",
    "y_2' = \\beta(t) y_1 y_3 - \\frac{y_2}{\\lambda} \\\\\n",
    "y_3' = \\frac{y_2}{\\lambda} - \\frac{y_3}{\\eta}\n",
    "$$\n",
    "Boundary conditions:\n",
    "$$\n",
    "y(0) = y(1)\n",
    "$$\n",
    "Parameters:\n",
    "$$\n",
    "\\mu = 0.02, \\lambda = 0.0279, \\eta = 0.01, \\beta(t) = 1575 (1 + \\cos 2 \\pi t)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mu = 0.02\n",
    "l = 0.0279\n",
    "eta = 0.01\n",
    "def fun_measles(x, y):\n",
    "    beta = 1575 * (1 + np.cos(2 * np.pi * x))\n",
    "    return np.vstack((\n",
    "        mu - beta * y[0] * y[2],\n",
    "        beta * y[0] * y[2] - y[1] / l,\n",
    "        y[1] / l - y[2] / eta\n",
    "    ))\n",
    "\n",
    "\n",
    "def bc_measles(ya, yb):\n",
    "    return ya - yb\n",
    "\n",
    "x_measles = np.linspace(0, 1, 5)\n",
    "y_measles = np.full((3, x_measles.shape[0]), 0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I use a custom tolerance to make the plot match the one from the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Iteration    Max residual    Total nodes    Nodes added  \n",
      "       1          1.39e-02           5              4       \n",
      "       2          2.41e-02           9              8       \n",
      "       3          1.45e-02          17             13       \n",
      "       4          1.95e-02          30             14       \n",
      "       5          9.14e-03          44             18       \n",
      "       6          1.69e-04          62              0       \n",
      "Solved in 6 iterations, number of nodes 62, maximum relative residual 1.69e-04.\n"
     ]
    }
   ],
   "source": [
    "res_measles = solve_bvp(fun_measles, bc_measles, x_measles, y_measles, verbose=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_measles_plot = np.linspace(0, 1, 100)\n",
    "y_measles_plot = res_measles.sol(x_measles_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1093bb7f0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEACAYAAACkvpHUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lGX28PHvCRCaSJMmJSCCCDYsqEiJAgKiAsoqoBRR\n4bcuu667IrLrCq6KAmthV31ZV1SwgUpvSo2K0hGVJiCKFOkgPZDkvH/cE4hhUiYzmWfK+VzXyExy\n38+cjMmcubuoKsYYY0x+JXgdgDHGmOhiicMYY0xALHEYY4wJiCUOY4wxAbHEYYwxJiCWOIwxxgQk\nJIlDRNqJyHoR2SAiA3Mo828R2Sgiq0TkirzqisjlIrJIRL4WkaUicnUoYjXGGBOcoBOHiCQArwBt\ngUZANxFpkK1Me6CuqtYD+gGj8lF3ODBYVRsDg4ERwcZqjDEmeKFocTQBNqrqFlU9BYwDOmYr0xEY\nC6CqS4CyIlIlj7oZQFnf/XLA9hDEaowxJkhFQ3CN6sDWLI+34RJCXmWq51H3EeBTEXkBEKBpCGI1\nxhgTJK8GxyUfZX4PPKyqtXBJ5M3CDckYY0x+hKLFsR2oleVxDc7uVtoO1PRTJjGXur1U9WEAVf1Y\nREb7e3IRsc22jDGmAFQ1Px/izxKKFscy4EIRSRKRRKArMDVbmalATwARuQ44qKq7cqg7xVdnu4i0\n9NVpBWzIKQBVtZsqgwcP9jyGSLnZa2Gvhb0Wud+CEXSLQ1XTRaQ/MBuXiEar6joR6ee+ra+r6kwR\nuUVENgFHgftyqbved+kHgX+LSBHgBNA32FiNMcYELxRdVajqJ8BF2b7232yP++e3ru/rXwG2dsMY\nYyKMrRyPIcnJyV6HEDHstTjDXosz7LUIDQm2r8trIqLR/jMYY0y4iQjq4eC4McaYOGKJwxhjTEAs\ncRhjjAmIJQ5jjDEBscRhjDEmIJY4jDHGBMQShzHGmIBY4jDGGBMQSxzGGGMCYonDGGNMQCxxGGOM\nCYglDmOMMQGxxGGMMSYgljiMMcYExBKHMcaYgFjiMMYYExBLHMYYYwISksQhIu1EZL2IbBCRgTmU\n+beIbBSRVSJyRX7qisgfRWSdiHwnIs+HItZgqMLmzfDGG9CjBzRvDhddBOXLQ5ky0KgR3HIL/PGP\nMH8+ZGR4HbExxoRe0EfHikgCsAFoBewAlgFdVXV9ljLtgf6q2kFErgVGqup1udUVkWTgb8Atqpom\nIuep6l4/z1/oR8f+8AOMHg3vvw8nT8KNN0JysksalSu7W0ICbN0KW7bA2rXw3ntw8CD07AkPPQTV\nqhVqiMYYE5Bgjo4tGoLnbwJsVNUtvmDGAR2B9VnKdATGAqjqEhEpKyJVgDq51P098LyqpvnqnZU0\nClNGBkyZAq+8At9951oY06bBJZeA5PBSlysHl14Kt94Kjz0Gq1a5hHPFFfDii9C9e851jTEmWoSi\nq6o6sDXL422+r+WnTG516wMtRGSxiCwQkatDEGueMjLg44/dm/0zz0Dfvq4l8cILLikE8sZ/xRXw\nn//AzJnw/PPQuTPs3Fl4sRtjTDiEosVREPl5+y0KlPd1aV0DfAhc4K/gkCFDTt9PTk4mOTm5QEHN\nmwd/+QskJsLQodChQ2haCFddBcuXwz//CVdf7cY/6tcP/rrGGJNfKSkppKSkhORaoUgc24FaWR7X\n8H0te5mafsok5lJ3GzARQFWXiUiGiFRU1X3ZA8iaOApi82Z49FHXtTRiBNxxR+i7lIoXh2efhbp1\noVUrWLAALrwwtM9hjDE5yf6h+qmnnirwtULRVbUMuFBEkkQkEegKTM1WZirQE0BErgMOququPOpO\nBm7y1akPFPOXNIJx8qRrBTRpAtdc4wa177yzcMch+vSBf/wDbrrJJSxjjIk2Qbc4VDVdRPoDs3GJ\naLSqrhORfu7b+rqqzhSRW0RkE3AUuC+3ur5Lvwm8KSLfAan4Ek+oLF4MDzwAderA119DzZp51wmV\nvn0hPd3Nzlq82GZcGWOiS9DTcb0W6HTc48fhb3+DcePg5Zfhrru8m+n05JMucXzyiZvOa4wx4RLM\ndNy4ertascINVP/yC6xeDXff7e302CefhMOHYeRI72IwxphAxUWLIy0NnnvOTY0dORK6dQtTcPmw\neTNcey3MnQuXX+51NMaYeOH1AsCI9vPPcM89blbTypVQo4bXEf3WBRecWRy4fDmULOl1RMYYk7uY\n7qqaMMGtm7j1Vpg9O/KSRqZ774XLLnOzrYwxJtLFZFfViRNuId+nn7r9pa691qPgArBrFzRs6Fod\ndep4HY0xJtbZ4HgWP/wATZvC7t2uayoakgZAlSrw8MPw9797HYkxxuQuphLHxIlw/fVw333w0UdQ\ntqzXEQXmL3+BlBTX6jDGmEgVE11VJ08qgwa5zQk/+sitAo9Wr78OH3zg9rOynXSNMYUl7ruqWrWC\nNWvcOo1oThrgtiTZuRNmzfI6EmOM8S8mEkfr1jBjBlSs6HUkwStaFIYNg4ED7QRBY0zhGDMmuPox\n0VUV7T9Ddqpw5ZVua/f27b2OxhgTK06ehEcegTlzYOPGOO+qijUi8Oc/w0sveR2JMSZW7Njhjrze\ntg2WLQvuWpY4IlTXru7I2tWrvY7EGBPtvvrKjf/ecgtMmhT8jFNLHBGqeHF46CHbANEYE5z//hc6\ndYL//Q+eeCI0O3HbGEcE270bLroINmyASpW8jsYYE01SU+GPf4SFC2HKFKhX77ffj/vpuLGqcmV3\nIuGoUV5HYoyJJjt3uoPi9uyBJUvOThrBssQR4f78Z3jtNffpwRhj8rJsmRvPaNvWbfRapkzon8MS\nR4S75BJo1Mj9AhhjTG7eeQc6dHBnDw0eXHgni1riiAL33w9jx3odhTEmUqWlwaOPwlNPwYIFbjC8\nMIUkcYhIOxFZLyIbRGRgDmX+LSIbRWSViFyR37oi8lcRyRCRCqGINRp17Oj6KXfs8DoSY0ykOXDA\ntTK++QaWLnU9FIUt6MQhIgnAK0BboBHQTUQaZCvTHqirqvWAfsCo/NQVkRpAG2BLsHFGs1Kl4I47\n3NkixhiTad06d3REw4Zuf7sKYfp4HYoWRxNgo6puUdVTwDigY7YyHYGxAKq6BCgrIlXyUfclYEAI\nYox6PXu6/WVidOaxMSZAM2dCy5YwaJDbZaJoGA8CD0XiqA5szfJ4m+9r+SmTY10RuR3YqqrfhSDG\nqNe8ORw+7Jqjxpj4pQojRsADD8Dkye78oXALY476jVwXnYhISeBvuG6qPOsMGTLk9P3k5GSSk5OD\niy4CJSRAjx5ukPyKK/Iub4yJPSdOQN++biuiJUugZs38101JSSElJSUkcQS9clxErgOGqGo73+PH\nAVXVYVnKjAIWqOp43+P1QEugjr+6wAxgLnAMlzBqANuBJqq6O9vzx+zK8ew2bIAWLdwmZeFslhpj\nvPfLL9C5MyQlwVtvubHPYHi9cnwZcKGIJIlIItAVmJqtzFSgJ5xONAdVdVdOdVV1tapWVdULVLUO\nrgurcfakEW/q14c6deDTT72OxBgTTitWuEHwW2+FceOCTxrBCjpxqGo60B+YDawBxqnqOhHpJyJ9\nfWVmAj+KyCbgv8BDudX19zTk0b0VL3r2dIt8jDHxYfx4aNcOXn7ZbVIYCUdK2yaHUWb3brfvzK5d\nUKKE19EYYwpLRoZb/f3OO26TwssvD+31ve6qMmFUubL7BZo71+tIjDGF5cgR6NIFUlLcor5QJ41g\nWeKIQp07u8NYjDGxZ8sWuOEGKF/efUCsXNnriM5miSMKde4MU6e6/WmMMbHjyy/h+uvd2ow33nAH\nukUiSxxRqHZtN3974UKvIzHGhMrbb7sPhaNHu+MUImEQPCe2GiBKZXZXxeBaR2PiSno6DBzoBsA/\n+wwuvtjriPJms6qi1Jo17uD5n36K7E8mxpic/fordO/uVoR/9FH4NikEm1UVlxo2dP2fK1d6HYkx\npiB++AGaNnUrwT/5JLxJI1jWVRWlRFx31cSJcNVVwV3r519/ZvG2xSzdvpTlO5ZzPO045UuUp1yJ\nciSVTaJTg05cW+NaEsQ+ZxgTCikp0LUrPPkkPPSQ19EEzrqqotjixdCnD6xdW7D6X/78Jc8tfI4l\n25dwQ80baFK9Cdecfw1lipfhwPEDHDxxkHV71/Hx2o85fPIwXS7uwl+b/pUa59YI7Q9iTBx5/XX4\nxz/c+TqtWnkXRzBdVZY4olhGBtSo4QbU6tXLf73F2xYzYM4AdhzewYCmA+h9RW9KFM19GfraPWsZ\ns2oMo78ezZ+u/ROPNn2UUsU83jDHmCiSlgZ//avba27atMD+ZguDJY4o/xmC0acPNG4Mf/xj3mVP\npZ/i6c+f5vUVrzOizQi6XdqNogmB9VZuObiFx+Y+xqKti3j1lle57aLbChi5MfHj4EG4+253f/x4\nKFfO23jABsfjWvv2bmAtLxv2beCGN29g2Y5lfN3va3pc3iPgpAGQVC6J8V3G807nd/jDzD8waO4g\n0jJsJaIxOdm4Ea67Dho0gBkzIiNpBMtaHFHuwAE3K2P37pw3PVzw4wLu/vhunmz5JH+45g9IiObv\n7jm6h24TuqEoH9z5AZVLR+DeCMZ4aN48N9326afdAUyRxFoccax8ebjsMvj8c//fH7d6HF0ndOXD\n331I/yb9Q5Y0ACqVrsSn937KtdWvpcn/mvD93u9Ddm1jot1rr8E997iuqUhLGsGy6bgxoF07mDUL\nbr75zNdUlRcXvcjIJSOZ13Mel1S+pFCeu0hCEYa2Gkr9ivW5ccyNzOg+g8bVGhfKcxkTDU6dcluG\nLFjg9p6qW9friELPWhwxwN84xz8/+ydvrnqTr+7/qtCSRla9r+jNK7e8Qtt32/LFli8K/fmMiUQH\nDri/xx9/hEWLYjNpgCWOmNC4Mezf77YfAXj282cZv2Y8C3otCOuaizsuvoP373yfOz+8k7mb7cAQ\nE1++/94d73rZZW66bdmyXkdUeGxwPEb07Om2Lzh06XBGfz2alF4pVCtTzZNYFv68kDvG38GEuybQ\nPKm5JzEYE05z5sC998LQoXD//V5Hkz+2jiPKf4ZQeP99eG7+SI5f9h8+6/0Z1c+t7mk8czfPpfuE\n7kzvPp0m1Zt4GosxhenVV92sqQ8/hBYtvI4m/zyfVSUi7URkvYhsEJGBOZT5t4hsFJFVInJFXnVF\nZLiIrPOVnyAi54Yi1lj1a9K7rCn7ArO6zfM8aQC0vqA1b3V8i9s+uI1VO1d5HY4xIXfqlNtn6rXX\n4KuvoitpBCvoxCEiCcArQFugEdBNRBpkK9MeqKuq9YB+wKh81J0NNFLVK4CNwKBgY41VszbOYsii\nv3LJqk/Y+l2S1+Gc1qF+B1695VVuee8WNu7b6HU4xoTM/v1uNuOWLW4Q/IILvI4ovELR4mgCbFTV\nLap6ChgHdMxWpiMwFkBVlwBlRaRKbnVVda6qZvjqLwZsZz0/Fm1dRM/JPZl892Q6N2vIp596HdFv\ndWnYhaeSn+Lmd29mx+EdXodjTNDWr3eD4I0buyOcz43DvpBQJI7qwNYsj7f5vpafMvmpC9AHmBV0\npDFm3Z51dB7fmbGdxnJ9zetp1Qrmz/c6qrM9eNWD9L2yL23fbcuB4we8DseYAps9G1q2hEGD4F//\ngiJFvI7IG14tAMz3gIyI/B04parv51RmyJAhp+8nJyeTHAfnqW4/tJ1277VjRJsRtK/XHnD74axf\n7+aSly/vcYDZPN7scXYf3c1tH9zG7B6zbWddE1VU4ZVX3Kypjz+G5lE4WTAlJYWUlJSQXCvoWVUi\nch0wRFXb+R4/DqiqDstSZhSwQFXH+x6vB1oCdXKrKyK9gQeBm1Q1NYfnj7tZVQdPHKT5W82599J7\nGdjst3MR2raF3/8eOnXyKLhcZGgGvSb34sDxA0y6exLFihTzOiRj8nTqlNt9+ssvXddUnTpeRxQa\nXs+qWgZcKCJJIpIIdAWmZiszFegJpxPNQVXdlVtdEWkHDABuzylpxKMTaSfoNK4TN9W+icdueOys\n77dq5TZWi0QJksCbt79Jhmbw4LQHibeEb6LPvn3uw9j27S5xxErSCFbQiUNV04H+uFlQa4BxqrpO\nRPqJSF9fmZnAjyKyCfgv8FBudX2X/g9wDjBHRFaKyGvBxhrt0jLS6D6hO1XOqcJL7V7yu2HhTTdF\n5jhHpmJFivHR7z7i+33fM3Cu35nbxkSEdevcIPhVV8HkyfE5CJ4TWwAYJVSVB6c9yNZDW5nWbRqJ\nRRL9lktPh0qVYM0aqObNwvF82X98P83fak6vy3v5bTkZ46VPPnG7MQwfDr17ex1N4fC6q8qEwaB5\ng1i9ezUT7pqQY9IAN8sjOTmyWx0AFUpWYPa9sxm1fBSjlo/yOhxjADcIPnIk3HcfTJwYu0kjWLat\nehQY/uVwpn4/lS/u+4JzEs/Js3xmd9U994QhuCBUP7c6c3rMoeXbLTkn8Rzuvexer0MycezkSejf\n3y3oW7QIatf2OqLIZS2OCPfy4pf574r/MrvHbCqWqpivOpkD5NHQg1e3Ql0+vfdTHp39KJPXT/Y6\nHBOn9u1z59ns3Om2D7GkkTtLHBHs1aWv8u8l/w54e/QGDdynpx9/LMTgQqhR5UbM6D6DftP7MfX7\n7BPyjClcmYPg114LkyZBmTJeRxT5LHFEqFHLRzHiqxHM7zWfWmVrBVRXxHVXReq0XH+uOv8qZnSf\nwYPTHrSWhwmbTz5xK8GfeAKGDYvfleCBssQRgUYuHslzC59jXs951C5Xu0DXiOT1HDm5+vyrmXXP\nLPpN78eEtRO8DsfEMFV4+WU3CD5pkg2CB8qm40YQVeXJBU/y4doPmdNjTsAtjay2bIEmTVyfrZ/l\nHhHt61++pv177RnRZgQ9Lu/hdTgmxmQOgi9e7FaCx+t4RjDTcW1WVYRIz0in/8z+LN2xlIX3LaRS\n6UpBXS8pCUqWdMdZNmiQR+H0dHfu7Lp1brOrHTvcvtH798PRo1CsGBQvDiVKuMUhtWq5W/360LAh\nFA3tr1Hjao2Z32s+7d9rzy9HfmFA0wF+FzsaE6h9++DOO91ivi+/tPGMgrIWRwQ4nHqYnpN7cvDE\nQaZ0ncK5xUOzRLVXL3ecbL9+2b5x4oSbb/j55+62ZAlUrAgXX+xuNWtChQruVrq0+4iWmgrHj8Mv\nv8DWra5Js3YtbNsGl17qmjetWsGNN4bsr3H7oe20f689N9W5iRfbvkiCWM+qKbh16+C221ziGDrU\nxjPs6Ngo/hk2H9jM7R/czvU1rueVW16heNHiIbv2m2/C3LnuWFn274fp02HKFPfFhg3dqGCLFi67\nlCtXsCc5dAi+/tq1++fMcUmocWO4/Xb3Fxrk5j4HTxyk07hOVChZgTGdxlCmuH1ENIHLXAk+YoT7\nQGUscURt4pi3eR73TLyHf7T4Bw9d81DIu2M2rz7GiOZTea3Ze8jnn7upVh07QocObl+SwnDsGHz2\nmdvcZ9Ik13q56y7o3t3dL4DUtFQe/uRhPtvyGZPunkSD8/LqezPGyVwJPmyY2w79hhu8jihyBJM4\nUNWovrkfIbqkpqXqoLmDtOq/quq8zfNCe/GMDNWFC1X79NGMcuV0QfG2+svwsaqHDoX2efLj1CnV\nefNUH3xQtUIF1ZtuUn3rLdXDhwt0uTdWvKGVhlfSiWsnhjZOE5NSU92v3qWXqv74o9fRRB7fe2eB\n3netxRFm6/as495J93J+mfN547Y3qHJOldBceNcuGDsWRo9206j69IEePbjnr1W58UZ44IHQPE2B\nnTjhusrGjIGFC+GOO1yMTZsGNO1r2fZl/O6j39Hmgja80PaFkI0Hmdiydy906eIGwd97zwbB/bFN\nDqPAibQTDP1iKC3ebkG/q/oxtevU4JNGejrMmuXGEho0cDOi3nzTDVoPGABVq9Kypes58lyJEu4v\nedo0F99FF8H997vB+OHD3bzhfLim+jV8+/tvUZTLR13Ogh8XFHLgJtqsXWsrwQtdQZsqkXIjCrqq\npn0/TeuOrKsdP+ioP+z/IfgLbtyo+re/qVavrnr11aqjRqn++qvfouvXq9as6XqwIk5GhuqXX6r2\n6aNarpzq7berTpjg+hjyYcaGGVr9heraZ3If/eXwL4UcrIkGs2apVqqk+vbbXkcS+Qiiq8rzN/5g\nb5GcOL7Y8oW2GdtG6/+nvs7aOCu4i+3f7xJEs2buL+ORR1S//TbPahkZqlWqqG7eHNzTF7rDh1Xf\nfFO1ZUvV885T7d9fddGiPDPeweMHdcDsAVpxWEUd+vlQPX7qeHjiNRElI0P1pZdUq1Z1Q3wmb8Ek\nDhvjCDFVZf6P83nmi2fYcnALg5oNotcVvXI9QyNHhw65cYGPPnL7pLdtCz16uH8T83+9u++G9u2j\naFuFH3904zUffOAOfO7a1f0Ql16a43jIpv2bGDBnAMt3LOeR6x7hwSsftKm7ceLkSfjDH9xM8Hhe\nCR4om44bAT/DnqN7GPPNGN5Y+QZFEoow8IaBdLukG8WKFAvsQjt2wMyZbiwgJQWaN3djA506FXit\nxWuvwbJl8NZbBaruHVW3RuSDD1zyLFoUOnd2t2uv9buCa8WOFQz7chgLflrA/131fzx41YNBbd1i\nIlvmSvCyZeHdd208IxCWODz6GfYc3cPU76cy+fvJLPx5IZ0adOKBxg/QtGbT/K/JOHrUzTJasABm\nz3Zbf7RtC7fe6m5lywYd55o1bsXs5s1BX8o7qrBqlRvtnDLFrVhv3RratXMr1mv9Njls3LeRlxa/\nxPg142lctTG9r+hNpwad8nUQlokOa9e63+suXWwleEF4njhEpB3wMm6W1mhVHeanzL+B9sBRoLeq\nrsqtroiUB8YDScBPwF2q+quf64YtcRw9eZSvtn7FZ1s+Y8FPC1i9ezVt67alc4POdKjfIe+poaqu\nG2bpUteuXrwYvvsOrrzSbdXRujVcf33I935ShcqVYcWKs95fo9f27S7RfvKJa5mVLOlaZ82awTXX\nuG6tYsU4kXaCqd9P5e1Vb/PFz19wQ80b6FCvA+3rtadu+bq2B1aUmjXLrQC3leAF52niEJEEYAPQ\nCtgBLAO6qur6LGXaA/1VtYOIXAuMVNXrcqsrIsOAfao6XEQGAuVV9XE/zx/yxKGq7Di8gw37NvDd\n7u/4eufXrPxlJZv2b6Jx1cYk106mZVJLmic1p0TREmdfIC3N7eW0caO7rV3rEsR337m2dJMmZ+YL\nNmkCpUqFNH5/7rzTLZ2I9ONkC0QVNmxw+259+SUsX+4S9CWXuATSqBE0bMjhOtWZnbqW6ZtnMfuH\n2aRnpNO0ZlOur3E9jas15pLKl1CldBVLJhFMfSvBhw93K8GbNvU6oujldeK4Dhisqu19jx/HjdYP\ny1JmFLBAVcf7Hq8DkoE6OdUVkfVAS1XdJSJVgRRVPWuviUASR4ZmcCj1EAeOH2Dvsb3sPrqbPcf2\n8MvhX9h6aKu7/bqVTfs3UaZ4GepVqEfDSg25stqVNK50GZeWqk2JQ8fcvk/79sHu3W7h3a5dbmxi\n61b4+We3EWC1alCvHlx4odsX6tJL3a1i/o5/DbWRI13++u9/PXn68DtyxHVtrV7t+upWr3ZJfO9e\nqFULveACjlQux0+lT7E68QDfsZuVaVvZVzqBirUuolLVutSucAFJ5ZI4v8z5VCldharnVOW8UueF\ndD8xk3+Zg+BLl7pB8KQkryOKbl5vq14d2Jrl8TagST7KVM+jbhVV3QWgqjtFpHJOAYz9842kp50i\nPf0U6WknSTuVSlraKdJPpZJ2KpX0k6mkpaWSkZpKKSlGaSnBuQklOFdKkCTFuUQTOVeLUyajKKXT\nS1Lq1GUUPX4Cju6CIz/AoffcHkznnHNm19jy5aFKFXerWhUuu8z1A9WsCdWru23II0jz5nGUNMD9\nv2rWzN2yOn4cfvoJ+fFHymzbxqXbtnHp9u1021UK3VOa9N07Ye93JBxbSmqp4hwtWYTDxYXDxTLY\nVDSNlUVOcaJYAhklEtESJSCxOBQvTkLx4khiIglFE0lILI4USyShaFESihZDihQjoWhRpEhR379F\nSJAikJBAQkIRJCEBkQQkIQEkARFBJAESfH/TkoD4/gU3sUyy7RT8m8cJZ94LhFzeFyTHB/mW6/VD\n6Nhx+OhjOLc4vHI/7JoJu4K4nhJ8L0X2a4j6/15uH2xVM7JdE9esyqyn6q7l+1dPfy/DTY1F3TUU\nNCODDE1HVcnISHdlMjLIyEg/fUvPSCMjPZ2MjLSgfnavzuMoyG9bjq/+kpkbT//xXV+lIjecX4eE\nYokUK55IscSSFEssQbHEEhQvcQ4JicXdGEJi4m9vJUu61c0lSriuo9Klz9zKlnVvRAnRu9D+8svd\nsMCePYW3v2FUKFnyzPbx2QhZ/iDS0yl56BAlDx7kvCNHXAvm8GH06FFSj/7KsUP7OHH4ACePH+Hk\n8aOkHT9C+slUMlJPkHHyJBlHTkL6MTLS0yAtDc1Ih/R0ND0dMjLcm4Nm+O5z5vHp++7XXdS9BUnm\nG0a2+2R5/Jsvnv2AnL4jBX4PDc/YYloaHDsId5Zwf4a8Hfw1Jct/g7pGDpc46/ri7/uOZusaPeuR\nZPv39HMLIL4Z6pL5icLX1er7N8F9CBFJYNmvx1lx8FiWMgUXisSxHcg65FrD97XsZWr6KZOYS92d\nIlIlS1fV7pwCeHXDtgKGHj+KFHH9wQsXutmsJg9FirhWZfnyv/myACV8N1P4bBA8dC4Hsm5Z9/+C\nSB6h+Ai9DLhQRJJEJBHoCkzNVmYq0BNOj4kc9HVD5VZ3KtDbd78XMCUEsca1Fi3c+LExkU4VXnrJ\nbWc2ebIljUgTdItDVdNFpD8wmzNTateJSD/3bX1dVWeKyC0isgk3Hfe+3Or6Lj0M+FBE+gBbgLuC\njTXetWgBf/qT11EYk7usK8EXLbJB8EhkCwDjSGqqm9S1Y4fbbtqYSLN3r5s6Xq6crQQvbLatusmX\n4sXh6qvhq6+8jsSYs2Vuh3799bYdeqSzxBFnmjeHL77wOgpjfmvmTEhOhsGD4fnno3oCY1yw/z1x\nxgbITSQMm7noAAAX3ElEQVTJHAR/4AE3CN6zp9cRmfywMY44c+SIW6+4d69bsmKMV2wluLdsjMPk\n2znnuB1Qli71OhITz/buhTZt3K49X35pSSPaWOKIQ82bW3eV8c6aNW4QvGlTNwh+ju10H3UsccSh\nFi1sgNx4Y+ZMd4LAkCHw3HM2CB6tbIwjDu3bB3XquE1+Q3z0hzF+qcLLL7utQ2w79Mjg9e64JspU\nrOj6lFetcus6jClMJ0/CQw+544ttJXhssIZinLJpuSYcMgfB9+yxQfBYYokjTtkAuSlsmYPgmSvB\nbRA8dtgYR5zavt2dPbVnjw1QmtCbORN694YXXoAePbyOxvhj6zhMwKpXd0dNrFuXd1lj8iv7SnBL\nGrHJBsfjWOY4R6NGXkdiYoENgscPa3HEMRsgN6GSOQi+d68NgscDSxxxLHOA3IaITDCyDoJPnGiD\n4PHAEkccu+AC9++PP3obh4leM2a47dCHDLHt0OOJ/W+OYyLWXWUKRhVefBEefBCmTLFB8HhjiSPO\nWeIwgTp50s2aGjPGDYLb9iHxJ6jEISLlRWS2iHwvIp+KSNkcyrUTkfUiskFEBuZVX0Rai8hyEflG\nRJaJyI3BxGlyZgsBTSAyB8H37bNB8HgWbIvjcWCuql4EzAcGZS8gIgnAK0BboBHQTUQa5FF/D3Cr\nql4O9AbeCTJOk4OGDeHgQbcg0JjcZJ4J3rSpDYLHu2ATR0dgjO/+GKCTnzJNgI2qukVVTwHjfPVy\nrK+q36jqTt/9NUAJESkWZKzGj4QEaNkSPvvM60hMJMs8E9y2QzcQfOKorKq7AHxv9JX9lKkObM3y\neJvvawBV8qovIl2Alb6kYwpBcjKkpHgdhYlEmYPgDzxgg+DmjDxXjovIHKBK1i8BCjzhp3iwKwJ+\nU19EGgHPAW1yqzRkyJDT95OTk0lOTg4yjPiSnAz/+Y/XUZhIYyvBY0tKSgopIfqEGNQmhyKyDkhW\n1V0iUhVYoKoXZytzHTBEVdv5Hj8OqKoOy62+iNQA5gG9VHVxLjHYJodBysiAypXhm2/cHlbG7NkD\nd94JFSrAu+/aeEYs8nKTw6m4wWuAXsAUP2WWAReKSJKIJAJdffVyrC8i5YDpwMDckoYJjYQENy3X\nxjkMwOrVbhD8hhtsENz4F2ziGAa0EZHvgVbA8wAiUk1EpgOoajrQH5gNrAHGqeq63OoDfwDqAk+K\nyNcislJEzgsyVpMLG+cw4FaC33QTPPWUDYKbnNl5HAaAb7+FLl1gwwavIzFeyBwEf+EFmDDB7Ttl\nYpudOW6CdsklblHX9u02zhFvUlPh97+HlSth8WKoVcvriEyks4aoAWw9R7zas8etBN+/HxYutKRh\n8scShznNxjniS+YgeLNmNghuAmOJw5xmiSN+zJzpBsH/+U8YOtQGwU1gbHDcnJaRAZUquYFyG+eI\nTZlngv/rXzYIHu+8XMdhYkjmOIe1OmJT5nboY8e6QXBLGqagLHGY32jVCubN8zoKE2p79kDr1jYI\nbkLDEof5jdatYe5cO4c8lmRdCT5hgg2Cm+DZOg7zG/Xru6SxcaO7b6LbjBlw331ucd+993odDdSu\nXZstW7Z4HUZcSUpK4qeffgrpNS1xmN8QOdPqsMQRvbKuBJ8yJXLGM7Zs2YJNZgkvkQKNf+fKuqrM\nWTITh4lOqalw//3wzjs2CG4Kh03HNWfZuRMuvtidL12kiNfRmEDs2QN33AHnnecSR6SNZ/imgHod\nRlzJ6TW36bgmpKpWhRo1YPlyryMxgcgcBG/RwgbBTeGyxGH8su6q6JK5HfrTT8Ozz9pKcFO47NfL\n+GWJIzqougHwvn1h6lS45x6vIzLxwMY4jF9Hjrguq127oHRpr6Mx/qSmujPBV6xwSSMaFvXZGAdM\nmTKFNWvWUKRIEc4//3x69OgRUFlVpXz58iQkJJx+LW+++WbGjx/v9xqFMcZhicPkqEUL+PvfoW1b\nryMx2UX6IHhO4j1xHDp0iBtvvJEVK1YAcP311zN9+nQqVqyY77KHDh1i0aJFNG3alISEBCZPnkyb\nNm24+OKL/T6nDY6bsLLuqshkg+DR6/PPP6dRo0anH19++eUsWLAgoLIlSpSgc+fO1K5dm3PPPZdi\nxYrlmDQKiyUOk6PWreHTT72OwmQ1fTrceKPbDt0GwaPPtm3bKFeu3OnH5cqVY+PGjQGVrVatGiVL\nlgRg1KhR9OnTp3CD9iOoXzsRKS8is0XkexH5VETK5lCunYisF5ENIjIwv/VFpJaIHBaRvwQTpymY\na6+FHTtg61avIzGqbiv0fv1g2rTI2D4kFm3bto2JEyfSrVs3AE6dOkWbNm1Cdv0DBw5QokSJ048T\nExM5cuRIgcoeOHCAffv2Ubx48ZDFl1/BbjnyODBXVYf7EsIg39dOE5EE4BWgFbADWCYiU1R1fT7q\nvwDMDDJGU0BFirjxjZkz3RuW8UbWM8EXLYqOQfBghGKHjIIOo6xfv54mTZowcuRIABYtWkTt2rWZ\nOHEi33//PYMGDTqrzvDhwzlx4kS251dEhF69epGUlHT662XKlGH//v2nHx8/fpyqVav6jSWvsuPH\njw97F1WmYBNHR6Cl7/4YIIVsiQNoAmxU1S0AIjLOV299bvVFpCOwGTgaZIwmCB06wLhxlji8kjkI\nXqmS2w49HsYzvBw7b926Nc8++yz3+OY1z5s3j5tvvpmrrrqK1atX+63z2GOP5fv6devWZXmWlbX7\n9u3jyiuvLFDZ+fPn07Nnz3w/dygF20NaWVV3AajqTqCynzLVgaydHdt8XwOokq1+FQAROQd4DHgK\nCP0OXSbf2rVzBztl+0BlwmD1amjSxB2u9fHH8ZE0IsGSJUto1qwZ4N6cW7duHbJrt2zZkpUrV55+\nvHLlSlq1agXA5s2bfzP7KbeyABs3bjw91hFuebY4RGQOvjf0zC8BCjzhp3iwnxUyfP8OBl5S1WO+\nnR1zTR5Dhgw5fT85OZnk5OQgwzCZKlSAyy5zyaNdO6+jiR/Tp0OfPu6YV1vUF16dOnVi+vTpzJ07\nl7S0NMqXL8+hQ4dCcu1SpUrx2GOP8cwzz6CqDBgwgMqV3eftLl26MHr0aBo3bpxnWYCKFStSPYAz\nnlNSUkgJ1fGeqlrgG7AO12oAqAqs81PmOuCTLI8fBwbmVh/4HNdNtRk4AOwFHsohBjWFa+hQ1f79\nvY4iPmRkqI4YoXr++aqLFnkdTehF+t/rvHnzdNCgQaqqOmTIEH3//fdVVfWnn37SIUOGeBlageX0\nmvu+XqD3/mC7qqYCvX33ewFT/JRZBlwoIkkikgh09dXLsb6qtlDVC1T1AuBlYKiqvhZkrKaAOnRw\neyHF8bqtsEhNda2M995z26Ffd53XEcWfihUrUq9ePd59910uuugiunXrxpEjR/j4449ZsWIFa9as\n8TrEiBDUynERqQB8CNQEtgB3qepBEakG/E9Vb/WVaweMxI2pjFbV53Orn+05BgOHVfXFHGLQYH4G\nkzdVN5Nn9my33boJvd274c473SD4O+/E7jYv8b5y3Au25YgfljjCo18/qFcPHn3U60hiz3ffwe23\nQ/fubnfbWF7UZ4kj/GzLEeOZDh3ceg4TWtOmue3Qn33WVoKb6GEtDpMvR4+63XK3bYOyfvcHMIHI\n3A79pZdg4kS3Sj8eWIsj/KzFYTxTurRbTzBtmteRRL/sg+DxkjRM7LDEYfLtd7+Djz7yOorotns3\ntGoFhw65leA1a3odkTGBs8Rh8q1jR1iwwL3pmcB9951rXSQnuwQcqzOnTOyzxGHyrVw5aN7cuqsK\nIusg+DPP2CC4iW7262sCYt1VgVGFESPg//7PbSPSvbvXERkTPJtVZQJy4AAkJbnZVeee63U0kS01\n1SWMVavcmeA2nmGzqrxgs6qM58qXh2bN3KdnkzMbBDexzBKHCdhdd1l3VW5sENzEOuuqMgHL7K7a\nvh3KlPE6msgybZpbo/Hyy7Yduj/WVRV+1lVlIoJ1V50t+yC4JY3oVNhnjseKYI+ONXGqWze3i6vv\n7yuuZR0EX7zYxjOCJU8Ff+inDi5Yq8bfmeMAkyZN4ttvv+W2227L8ajXeGJdVaZAjh2DGjXgm2/i\n+41y9253JnjlyrG9HXqoRENX1bPPPkulSpXo27cvgwcP5pxzzqFly5ZcfPHF9OvXj/fff9/rEANi\nXVUmYpQqBV27wttvex2Jd7791g2C33ijOxPckkZsyH7m+AMPPECTJk3Ytm0bderU8Ti6yGAtDlNg\nK1e6T9ubN8ffSujMQfCRI21RXyCiocXx5ptvsnfvXkqUKMEHH3xwurvqueee4+GHH6ZUqVIeRxgY\na3GYiHLllVChAsyb53Uk4aMKw4fbSvBYNX/+fDZt2sRjjz3GgQMH+NOf/gTAtGnT6N+/P9u3b/c4\nwshgLQ4TlNdeg88+g/HjvY6k8KWmQt++sHo1TJnixnhMYCK9xfHNN9+wcuVKihUrRtGiRenatSuT\nJk3iueeeo1y5crRs2ZK///3vXocZkIg7OlZEygPjgSTgJ9yZ4b/6KdcOeJkzZ44Py6u+iFwGjALO\nBdKBa1T1pJ9rW+Lw0MGDULs2bNoE553ndTSFZ9cu1y13/vluXMfGMwom0hNHLIrErqrHgbmqehEw\nHxjkJ7gE4BWgLdAI6CYiDXKrLyJFgHeAvqp6CZAMnAoyVlMIypVz52W/+67XkRSeVaugSRNo08a1\nrCxpmHgXbOLoCIzx3R8DdPJTpgmwUVW3qOopYJyvXm71bwa+UdXVAKp6wJoVkev+++H11yEjw+tI\nQm/iRJcwRoyAIUPibxKAMf4E+2dQWVV3AajqTqCynzLVga1ZHm/zfQ2gSg716wOIyCcislxEBgQZ\npylELVpAiRIwY4bXkYSOKjz9NDz8MMya5fbnMsY4ea4cF5E5QJWsXwIUeMJP8WBbBZn1iwI3AFcD\nJ4B5IrJcVRcEeX1TCETgb39zhxTdeqt7HM2OHYPeveHnn2HpUqhWzeuIjIkseSYOVc1xoxYR2SUi\nVVR1l4hUBXb7KbYdqJXlcQ3f1wB25lB/G/C5qh7wPc9M4ErAb+IYMmTI6fvJyckkJyfn9WOZEOvc\nGZ54wh0te9NNXkdTcFu3uiNyL7kEUlJcS8qYWJCSkkJKSkpIrhXsrKphwH5VHSYiA4Hyqvp4tjJF\ngO+BVsAvwFKgm6quy6m+iJQD5gLNgDRgFvCiqs7yE4MNf0SIMWNg7NjoXdexcKHrkvrLX+Cvf43+\nllMksllV4ReJ03ErAB8CNYEtuOm0B0WkGvA/Vb3VV64dMJIz03Gfz62+73vdgb8BGcAMVT1rxpav\nnCWOCHHqFFx4oZt5dN11XkcTmP/9z7WYxoyBdu28jiZ2WeIIv4hLHJHAEkdkefVVmD3bLZCLBidP\nwiOPuFbS1KlQv77XEcU2SxzhZ4nDD0sckeX4cbjgAjcT6YorvI4mdzt3wu9+584XeecdKFvW64hi\nnyWO8IvEBYDG/EbJkjB4MPzhD5G9rmPJErjmGmjdGiZPtqRhTCAscZiQ69sX0tLgrbe8juRsqjBq\nFNx2m+tWGzzYFvUZEyjrqjKF4uuv3SDzmjWRs4fV0aNnTuqbMMHGM7xgXVXhZ11VJmo0buyOlR04\n0OtInA0b3EwvEXe8qyUN44+dOZ4/ljhMofnnP+HTT+HLL72NY+xYuOEG6N/fTbe1TQojnEjwtwLK\nPHN8x44dwJkzxz/++GOefvppVq5cGZIfMdpZV5UpVBMnwoABbjA63F1Whw/DQw/BihVubcmll4b3\n+c3ZoqGrKvuZ4xdeeCE1atRg7969pKen07VrV69DDIh1VZmoc8cdbjX2HXe4g5DCZeFC111WsiQs\nX25Jw+Rf9jPHb731Vs4//3yWLl3KnXfe6XF0kcFaHKbQZWRAly5w7rluplVhbuVx/LhbAf7BB+50\nwk7+Nvo3nomGFkdOZ44vWrSIadOmMXToUI8jDIy1OExUSkhwC+y+/dad111YPvvMtTJ27HDPZUnD\nBMrfmeMDBw5k3bp1lCxZkg0bNngdYkSwFocJm+3b3SD1vffCU09BkSKhue7PP7txlMWL4aWXXLeY\niUyR3uLwd+b44sWL2b17N2vXruW2226jUaNGXocZENtyxA9LHNFl9264+24oXhzefx8qVCj4tfbt\ng5EjXZdU//7w2GNQqlToYjWhF+mJIxZZV5WJepUrw5w57ryLq6929wN9H9m61W1MWK+e65Zavtwd\n62pJw5jwsBaH8cykSfCPf7jB8ocfhnvucbOg/Nm2ze0pNWkSrFzpzjl/5BGoXt1/eROZrMURftZV\n5YcljuimCvPnuy6nuXOhVi23u25SEvz6qxu/+Plnt11Ihw7upMG2ba11Ea0scYSfJQ4/LHHEjmPH\n4KefYPNm92+5ci6R1KoFNWpA0TwPOjaRzhJH+Fni8MMShzHRwxJH+NnguDHGGM9Z498YEzZJSUlI\nYW4dYM6SlJQU8msG1VUlIuWB8UAS8BNwl6r+6qdcO+BlXAtntKoOy62+iBQF3gCuBIoA76jq8znE\nYF1VxhgTIC+7qh4H5qrqRcB8YJCf4BKAV4C2QCOgm4g0yKP+74BEVb0MuBroJyK1gow15qWkpHgd\nQsSw1+IMey3OsNciNIJNHB2BMb77YwB/uwM1ATaq6hZVPQWM89XLrb4CpUWkCFAKSAUOBRlrzLM/\nijPstTjDXosz7LUIjWATR2VV3QWgqjuByn7KVAe2Znm8zfc1gCrZ6lfxff1j4BjwC64L61+qejDI\nWI0xxoRAnoPjIjKHM2/oAIJrETzhp3iwgw0Zvn+vBdKAqkBF4AsRmauqPwV5fWOMMcFS1QLfgHW4\nVgO4N/l1fspcB3yS5fHjwMDc6uPGRO7JUmc00CWHGNRudrOb3ewW+K2g7/3BTsedCvQGhgG9gCl+\nyiwDLhSRJFzXU1egm5/6vbPU/xm4CXhPRErjks9L/gIo6KwAY4wxBRPsdNwKwIdATWALbjrtQRGp\nBvxPVW/1lWsHjOTMdNzn86hfGngLaOh7qjdV9cUCB2qMMSZkon7LEWOMMeEVNVuOiEg7EVkvIhtE\nZGAOZf4tIhtFZJWIXBHuGMMlr9dCRLqLyDe+20IRudSLOMMhP78XvnLXiMgpEYnZ8wHz+TeSLCJf\ni8hqEVkQ7hjDJR9/I+eKyFTfe8V3ItLbgzALnYiMFpFdIvJtLmUCf98MZnA8XDdcgtuEW2FeDFgF\nNMhWpj0ww3f/WmCx13F7+FpcB5T13W8Xz69FlnLzgOnAHV7H7eHvRVlgDVDd9/g8r+P28LUYBDyX\n+ToA+4CiXsdeCK9FM+AK4Nscvl+g981oaXHktogwU0dgLICqLgHKikgVYk+er4WqLtYzW78s5sy6\nmViTn98LgD/i1gbtDmdwYZaf16I7MEFVtwOo6t4wxxgu+XktFCjju18G2KeqaWGMMSxUdSFwIJci\nBXrfjJbEkdsiwpzKbPdTJhbk57XI6gFgVqFG5J08XwsROR/opKr/D7cGKVbl5/eiPlBBRBaIyDIR\n6RG26MIrP6/FK0BDEdkBfAM8HKbYIk2B3jdtd9wYJiI3Avfhmqvx6mUgax93LCePvBTFbRx6E1Aa\nWCQii1R1k7dheaIt8LWq3iQidYE5InKZqh7xOrBoEC2JYzuQdZPDGr6vZS9TM48ysSA/rwUichnw\nOtBOVXNrqkaz/LwWVwPjxO3lfR7QXkROqerUMMUYLvl5LbYBe1X1BHBCRD4HLseNB8SS/LwW9wHP\nAajqDyLyI9AAWB6WCCNHgd43o6Wr6vQiQhFJxC0izP6HPxXoCSAi1wEH1bcPVozJ87Xw7SQ8Aeih\nqj94EGO45PlaqOoFvlsd3DjHQzGYNCB/fyNTgGYiUkRESuEGQ9eFOc5wyM9rsQVoDeDr068PbA5r\nlOEj5NzSLtD7ZlS0OFQ1XUT6A7M5s4hwnYj0c9/W11V1pojcIiKbgKO4TxQxJz+vBfAPoALwmu+T\n9ilVbeJd1IUjn6/Fb6qEPcgwyeffyHoR+RT4FkgHXlfVtR6GXSjy+XvxDPB2lmmqj6nqfo9CLjQi\n8j6QDFQUkZ+BwUAiQb5v2gJAY4wxAYmWripjjDERwhKHMcaYgFjiMMYYExBLHMYYYwJiicMYY0xA\nLHEYY4wJiCUOY4wxAbHEYYwxJiD/H3QIDJBTH800AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109373a20>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x_measles_plot, y_measles_plot[0] - 0.07, label='$y_1 - 0.07$')\n",
    "plt.plot(x_measles_plot, y_measles_plot[1], label='$y_2$')\n",
    "plt.plot(x_measles_plot, y_measles_plot[2], label='$y_3$')\n",
    "plt.legend(loc='lower right')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Flow in a vertical channel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equations:\n",
    "$$\n",
    "f''' - R[(f')^2 - f f''] + RA = 0 \\\\\n",
    "h'' + R f h' + 1 = 0 \\\\\n",
    "\\theta'' + P f \\theta' = 0\n",
    "$$\n",
    "Bondary conditions:\n",
    "$$\n",
    "f(0) = f'(0) = 0, \\: f(1) = 1, \\: f'(1) = 0 \\\\\n",
    "h(0) = h(1) = 0, \\: \\theta(0) = 0, \\: \\theta(1) = 1\n",
    "$$\n",
    "Parameters:\n",
    "$$\n",
    "R = 100, P = 70, A \\text { is unknown}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fun_flow(x, y, p):\n",
    "    A = p[0]\n",
    "    return np.vstack((\n",
    "        y[1], y[2], 100 * (y[1]**2 - y[0]*y[2] - A),\n",
    "        y[4], -100 * y[0] * y[4] - 1, y[6], -70 * y[0] * y[6]\n",
    "    ))\n",
    "\n",
    "\n",
    "def bc_flow(ya, yb, p):\n",
    "    A = p[0]\n",
    "    return np.array([ya[0], ya[1], yb[0] - 1, yb[1], ya[3], yb[3], ya[5], yb[5] - 1])\n",
    "\n",
    "x_flow = np.linspace(0, 1, 10)\n",
    "y_flow = np.ones((7, x_flow.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Iteration    Max residual    Total nodes    Nodes added  \n",
      "       1          1.30e+00          10             18       \n",
      "       2          2.42e-02          28             12       \n",
      "       3          3.47e-03          40              2       \n",
      "       4          7.51e-04          42              0       \n",
      "Solved in 4 iterations, number of nodes 42, maximum relative residual 7.51e-04.\n"
     ]
    }
   ],
   "source": [
    "res_flow = solve_bvp(fun_flow, bc_flow, x_flow, y_flow, p=[1], verbose=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found A=2.7606 --- matches with the value from the paper.\n"
     ]
    }
   ],
   "source": [
    "print(\"Found A={:.4f} --- matches with the value from the paper.\".format(res_flow.p[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_flow_plot = np.linspace(0, 1, 100)\n",
    "y_flow_plot = res_flow.sol(x_flow_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1094dd0b8>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4lNW1x/HvCheLXAuICkjQilqsIKjcRBjES0AFRBHE\n0kqrUCw9reep1drjIVqVaq3SFm9Yq7WKoKKCgOAFAqKCWKgoFwE5IARFRSteAAns88eOEkMuk2Rm\n9rwzv8/zzONM8uadxWuysrPftdc25xwiIpKZckIHICIiyaMkLyKSwZTkRUQymJK8iEgGU5IXEclg\nSvIiIhms0iRvZveb2TYzW1HBMTEzW25mb5nZ/MSGKCIi1WWV1cmbWU/gc+Ah51yHMj7fGHgFOMs5\nV2hmzZ1zHyUlWhERqZJKR/LOuUXAJxUcMhyY5pwrLD5eCV5EJE0kYk7+GKCpmc03s6VmNiIB5xQR\nkQSonaBzdAZOB+oDr5rZq8659Qk4t4iI1EAikvwW4CPn3C5gl5ktBDoCByR5M1OjHBGRanDOWXW+\nLt4kb8WPskwH/mpmtYCDgK7A7eWdKJMaon30EUyeDE8+Cf/6F7RqBV26wJFHQosW/lG/PuzdC0VF\nsGsXbNsGW7fCnDn5HHxwPqtWwcEHwwknQLdu/tG1KzRvHvpflzr5+fnk5+eHDiMt6Frsp2uxn1m1\n8jsQR5I3s8lADGhmZu8C44C6gHPOTXLOrTGzucAKYC8wyTm3qtoRRUBBAfz1r/Dii3DuuXDVVdCj\nB3z3u/Gf4+CDIT8fnIPCQvj3v2HJEpgwAV57DXJzoW9f/+jdGxo2TNa/RkQyWaVJ3jk3PI5jbgNu\nS0hEaeydd+DXv4Y33oBrroG//x0aN67ZOc2gdWv/OPdc/7GiIv+XwYsvwh13wPDh0L27//x550Hb\ntjX+p4hIltCK1zh89RX89rd+GqVrV1i1CkaNqlmCj8Vi5X6udm3/Ptde6xN9YSGMHg3LlsEpp8DJ\nJ8Mtt/hfOpmgomuRbXQt9tO1SIxKF0Ml9M3MXNTm5D/4AC68EJo0gXvvhcMPDxtPUREsXAhPPOHv\nBRxxBFxyCQwbBocdFjY2EUkOM6v2jVcl+QosXw7nnw8/+pGfP89Js797iopg3jx45BGYPt3ftB05\nEgYOhO98J3R0IsnRtm1bNm3aFDqMpMjNzWXjxo0HfFxJPgnmzPHJ/c47YciQ0NFU7ssv4amn4IEH\n/E3ciy+Gyy+HDgc0ohCJtuKEFzqMpCjv36Ykn2BLlvgbnNOn+xueUbNxIzz4IPztb346Z/RoGDoU\n6tULHZlIzSnJV/GcSvLf9vbbEIv5BHnOOaGjqZmiIpg9G+65B5Yu9VM5Y8b4On6RqFKSr5o0m2UO\na+tWyMuDm2+OfoIHX6UzYIBP9K++6hdlnXyyn7OfP9/X6ItIZlOSL7ZnDwwaBJdd5ke8meboo+FP\nf4J334X+/eGKK6BzZ3joIV8iKiKpMWvWLHbu3Jmy99N0TbEbbvCj3dmz/QKlTLdvn7+5fPvtsGYN\n/OpXvva/UaPQkYlULOrTNe3atWPZsmU0LGMZezKmaxLRoCzy/vUvmDjRl0xmQ4IHXw7av79/LFsG\nf/yjn6sfNQquvNL33RGRmpswYQIffvghubm59O7dm169epWZ4JMl66drdu6EESPgz3/2DcayUefO\n8Oij8Prr8OmncNxx8F//BZs3h45MJNp27NjB1KlTGThwID179mTevHmMHj06pTFkfZL/n/+BH/zA\nrxjNdkceCXfdBStX+sVUJ57oyy8zdN2JSNItWbKETp060aVLF9q3b8/nn39Oly5dUhpDVif5Zct8\nq+C77sqeaZp4HH443HqrLydt3tyP9EeNUrIXqYolS5YwYcIEioqKePrppwG46qqrUh5HVif5a6+F\n667Lrt7tVdG8Odx0E6xd6+foO3f2VTmFhaEjE6mYWWIeNdG1a1fq1avHL3/5SwYNGpSYf1g1ZG2S\nnz8f1q/3S/+lYs2awY03+pF9gwZ+g5Mrr4QPPwwdmUjZnEvMo6ZWr15N+/bta36iGsjKJO+c7wf/\n+99DnTqho4mO5s39NM7KlX5dwXHHwfXXw2efhY5MJP188MEHHHLIITXa1SkRsjLJP/20XwA0dGjo\nSKLp8MN9yenSpf6voXbt/E5ZWlQlst+SJUvo0aNH6DCyL8kXFcHvfgfjx6df6+CoOeoo+Oc/Ye5c\nmDULjj/e97mP8DoVkRpbtmwZY8aMYfHixQxNg5Fk1qW5Rx7xNxHPPjt0JJmjY0e/evauu/yN2lNP\n9Z08RbJRTk4OrVu3plmzZnTs2DF0OJW3NTCz+4FzgW3OuXK7k5vZKcArwFDn3JPlHBO0rYFzvvb7\nj3+Es84KFkZG27vXj+5/9zvfzXP8eGjTJnRUkkmi3tagIqG6UD4AVDjuNbMc4A/A3OoEkSqLFsHu\n3XDGGaEjyVy1asGll/pKnKOPhk6dYNw4v6mJiKRepUneObcI+KSSw34BPAF8kIigkmXiRPj5zzUX\nnwoNGvjKm+XLfcL//vdh6lTN14ukWo3TnZm1BAY55+4G0nbdaGEhPP88/PjHoSPJLm3awJQpfgpn\n/Hjo0wfeeit0VCLZIxFj2gnA1SVep2WinzTJ73uqVrph9Orlu30OGeIT/ZVX+mZoIpJciWg1fDIw\nxXzFf3Ogn5ntcc7NKOvg/Pz8b57HYjFisVgCQqjYV1/5JP/ii0l/K6lArVp+uuyii+C3v4X27f1G\nJkOHqneQSEkFBQUUFBQk5FxxbRpiZm2BZ5xzJ1Ry3APFx6VVdc2jj/o9W5Xk08srr/g9Z1u0gDvv\nhGOOCR2RRIGqa6qm0ukaM5uML408xszeNbORZjbazEaVcXhaXvk774SxY0NHIaX16OGncPr188+v\nv95XP4lI4mT89n//93/QpYvfpFt9atLX5s3+F/HatX5q7bTTQkck6apt27ZsytC+17m5uWzcuPGA\nj9dkJJ/xSf6WW2DDBrj33pS+rVSDc/DUU35Xqv79fTO0Jk1CRyUSXrIXQ0Xa1KlqRBYVZjB4sO9y\nmZPjd+yaPj10VCLRltEj+XXr/J/9hYW+skOiZcEC3++/Uyff5VKbi0u20ki+HFOnwoUXKsFHVe/e\n8MYbfu/ZDh20YlakOjJ6JH/CCb4zom7iRd9rr8HIkXDssXD33XDooaEjEkkdjeTLsGoVfPKJb3sr\n0deli994/bjjfGvjxx4LHZFINGTsSH7cONixA+64IyVvJym0ZInvQdSxo18DoY3YJdNpJF+Kc37+\ndtiw0JFIMnTt6rtbtm7tE/3s2aEjEklfGTmSX7ECBgzwC6HUEyWzLVjgR/V5eXDbbb7FsUim0Ui+\nlFmzfJJXgs98vXv7X+q7d/tSS207KPJtGZnk587VHq7ZpFEjeOABv7p54EC44Qa/YbuIZOB0zWef\nQcuW8P77UL9+Ut9K0tDWrX77wc8/h4cfhqOOCh2RSM1puqaEefP8jTkl+OzUsiXMmeM3J+nWDSZP\nDh2RSFgZl+Q1VSM5OX7nqblz/dTNj3/s/8ITyUYZleSd86O4vLzQkUg66NTJ96s/6CDo3Nk/F8k2\nGZXk16/3VRY/+EHoSCRd1K/v+9PfdJPfnOTPf1b/G8kuGZXk586Fs85S6aQc6KKLYPFieOQRX4Gz\nfXvoiERSI+OSvKZqpDxHHQWLFvm9ZE86ySd9kUyXMSWUu3fDIYf4Va7NmiXlLSSDzJjhe9X/5jfw\n3/+tv/4kvSV7I+/7zWybma0o5/PDzeyN4sciMzuhOoHU1MsvQ/v2SvASnwED/OrYxx6DQYN8x1KR\nTBTPdM0DQEVFiRuAXs65jsCNwH2JCKyqVDopVdW2Lbz0kv/vSSf5VsYimabSJO+cWwSUO85xzi12\nzn1a/HIx0CpBsVXJggXQp0+Id5Yoq1vXV9zccosfJEyapOobySyJvvF6GfBsgs9ZqZ074c03/cYS\nItUxZIi/KfuXv8BPfuK/p0QyQe1EncjM+gAjgZ4VHZefn//N81gsRiwWq/F7L13qa+MPPrjGp5Is\nduyxfp7+8suhRw+YNk29bySMgoICCgoKEnKuuKprzCwXeMY516Gcz3cApgF5zrl3KjhPUqprxo+H\nDz+E229P+KklCzkHEyfCjTfCgw/6RVQiIaWiQZkVP8p68zb4BD+iogSfTIsWaS9XSRwz+MUv4Mkn\n4bLLfLLfty90VCLVU+lI3swmAzGgGbANGAfUBZxzbpKZ3QcMBjbhfxHscc6VOTuejJH8vn1+j89V\nq+CwwxJ6ahG2boULL4QWLeChh3zvepFUq8lIPvKLoVau9DXP7wT5G0KywVdfwS9/CQUF8PTTfu5e\nJJWyup/8yy9rqkaSq25duPtuvzL2tNP89pIiUaEkLxKnyy+Hp56CUaPg5ptVTy/REPkkv2gR9Kyw\naFMkcU49FV57zU/bXHwxfPll6IhEKhbpJP/++77nyPe/HzoSySatWvkV1nXq+OmbzZtDRyRSvkgn\n+Zdfhu7d/XZvIqlUr56vthk2zO8l++qroSMSKVuk06Pm4yUkM7jqKt/vZuBAePjh0BGJHEhJXqSG\nzjkH5s+H//1fuPZaLZyS9BLZOvldu6BpU/joI/WskfTw4YdwwQV+T4OHH/b7y4okQlbWyb/5JrRr\npwQv6eOQQ+D556FxY+jVCwoLQ0ckEuEkv2wZdO4cOgqRbzvoIHjgAd+6uFs3bUQi4SnJiySYGVxz\nDUyY4DcimT49dESSzSKb5JcvV5KX9HbBBTB7Nlxxhd99SitkJYRI3njds8fPe37wATRokIDARJJo\n0yZfgdOnD9xxB9RO2FY9ki2y7sbr6tXQpo0SvERDbq4v912zBs4/H774InREkk0imeQ1VSNR07ix\nn7o55BDo3du35BBJhUgm+WXLoFOn0FGIVE2dOnD//X51bPfufqMbkWSLbJLXSF6iyAyuuw5uuMHP\n0S9YEDoiyXSRu/G6bx80aQIbN/oVryJR9eKLvl3xxIlw0UWho5F0VpMbr5G7z79+vV82rgQvUde3\nr18he+65fnXslVeGjkgyUaXTNWZ2v5ltM7MVFRzzFzNbZ2b/NrMTExvit2k+XjJJx46+8uZvf/Pb\nC6q5mSRaPHPyDwBnl/dJM+sHfM851w4YDdyToNjKpMoayTRt2sBLL8HSpXDJJbB7d+iIJJNUmuSd\nc4uATyo4ZCDwUPGxS4DGZnZoYsI7kG66SiZq2hSee84n+P79YceO0BFJpkhEdU0roOQGaIXFH0s4\n5zRdI5mrXj14/HE49ljV0kvipPzGa35+/jfPY7EYsVgs7q/dvBnq1oXDD098XCLpoFYtuPNOuPFG\nvyHO3Llw9NGho5JUKygooKCgICHniquE0sxygWeccx3K+Nw9wHzn3NTi12uA3s65bWUcW6MSyqef\nhvvug1mzqn0Kkci49164/nqYOVNTlNkuFb1rrPhRlhnAj4oD6Qb8p6wEnwhvvgknnJCMM4ukn9Gj\nfQ19Xp7fXlCkOiqdrjGzyUAMaGZm7wLjgLqAc85Ncs7NNrP+ZrYe+AIYmaxgV63yN6VEssXgwfDd\n78LQoX5kf/75oSOSqInUiteOHf2uO/rTVbLNsmW+XfFNN8FPfhI6Gkm1mkzXRCbJFxVBw4awfbv2\ndZXstHYtnHUW/PzncNVVoaORVMqKtgYbNviqGiV4yVbHHAOLFvlE//HHcPPNvuGZSEUi04Vy1Spo\n3z50FCJhtW4NCxfCCy/AmDGwd2/oiCTdRSrJH3986ChEwmveHObN89M3l1wCX30VOiJJZ5FJ8itX\naiQv8rWGDf1OU19+6Stwdu4MHZGkq8gkeU3XiHzbd74D06ZBo0bQrx989lnoiCQdRSLJ790Lb78N\n3/9+6EhE0kudOvDPf/p+N337+huyIiVFIslv3AgtWkCDBqEjEUk/tWrBPff4pmaxGGxLynpziapI\nlFBqPl6kYmZw661+INSrl99asHXr0FFJOohEktd8vEjlzGDcOKhfH047zSf6o44KHZWEFpkkX4WO\nxCJZ7de/9om+d29fT3/ssaEjkpAiMSevGnmRqhkzBn7/e+jTx3dvleyV9iP5fftg9WpV1ohU1aWX\n+jLLM8/0NfVq7Jed0j7Jb9rk979s1Ch0JCLRM2yYT/T9+sH06dCtW+iIJNXSPsnrpqtIzQwa5LfN\nHDAAnnwSevYMHZGkUtrPyWs+XqTm+veHRx7xm45ol6nsEokkr5G8SM2deSY8/jhcdBE8/3zoaCRV\nIpHkddNVJDFiMXjqKd+9cs6c0NFIKqR9kl+/3m+WICKJ0bOnvwn7ox/BrFmho5FkiyvJm1mema0x\ns7VmdnUZn29kZjPM7N9m9qaZXZqI4D7+2G/717x5Is4mIl/r3h2eecbvFztjRuhoJJkqra4xsxxg\nItAX2AosNbPpzrk1JQ77ObDSOTfAzJoDb5vZw865opoEt349HH20tjgTSYauXf1I/pxzwDkYODB0\nRJIM8ZRQdgHWOec2AZjZFGAgUDLJO6Bh8fOGwPaaJnjwSb5du5qeRUTKc/LJfqFU//5+4eH554eO\nSBItniTfCthc4vUWfOIvaSIww8y2Ag2AoYkIbt06P5IXkeQ56SR49lm/YMo5v9OUZI5ELYY6G1ju\nnDvdzL4HPG9mHZxzn5c+MD8//5vnsViMWAWdx9av9xshiEhyde7sq2369fOvlejDKigooKCgICHn\nMudcxQeYdQPynXN5xa+vAZxz7pYSx8wExjvnXi5+/SJwtXPu9VLncpW9X0ndu8Ntt8Gpp8b9JSJS\nA8uXQ16e34REUzfpw8xwzlXr7mQ81TVLgaPNLNfM6gLDgNL34zcBZxQHcyhwDLChOgGVpOkakdTq\n1MlP3fzsZ76eXqKv0uka59xeMxsLPIf/pXC/c261mY32n3aTgBuBB81sRfGX/cY5V6PdJj/5BL76\nym/7JyKp07nz/jl6M9/7RqIrrjl559wc4NhSH7u3xPP38PPyCaPySZFwOnfeX3VTqxacd17oiKS6\n0nbF69dJXkTCOOkkmDkTfvpTn/AlmtI6yatGXiSsU07xK2MvvVS9bqIqbZO8brqKpIeuXff3unnh\nhdDRSFWlbZLXdI1I+ujeHaZNg+HDIUHl25IiaZ3kNV0jkj5OOw0eewyGDIGXXgodjcQrLZP8p5/C\nzp1w6KGhIxGRkmIxePRRuOACePXV0NFIPNIyyat8UiR9nXEG/OMfvmvl669XfryEldZJXkTSU79+\ncN99cO658MYboaORiiSqQVlCqbJGJP0NHOhXpefl+aqb448PHZGUJS2T/Pr10KtX6ChEpDJDhvhE\nf9ZZvupGxRLpR9M1IlIjl1wCN9zg24Jv3Bg6GiktLUfymq4RiZaf/tRXxPXtCwsXQqtWoSOSr6Vd\nkt+xA774Ag4/PHQkIlIVY8fuT/QLFqgEOl2kXZLfsAGOPFLlkyJRdNVVfpB25pl+jr5p09ARSdrN\nyW/c6JO8iETTuHFw9tm+6mbHjtDRSNol+U2boG3b0FGISHWZwa23wskn+zr6L74IHVF2S7skv3Ej\n5OaGjkJEasIMJk70f5UPHgy7d4eOKHulXZLXSF4kM+TkwP33Q8OGcPHFUFQUOqLslHZJXiN5kcxR\nuzZMnuyrbkaOhH37QkeUfeJK8maWZ2ZrzGytmV1dzjExM1tuZm+Z2fzqBrRxo0byIpmkbl3fi/7d\nd+GKK8C50BFlF3OVXHEzywHWAn2BrcBSYJhzbk2JYxoDrwBnOecKzay5c+6jMs7lKnq/HTugZUv4\n7DOVUIpkmh07fAfLPn3gD3/Qz3hVmBnOuWpdsXhG8l2Adc65Tc65PcAUYGCpY4YD05xzhQBlJfh4\nbNrkp2r0P18k8zRqBM8+C7NmwfjxoaPJHvEk+VbA5hKvtxR/rKRjgKZmNt/MlprZiOoEo6kakczW\nrBk895y/ITtxYuhoskOiVrzWBjoDpwP1gVfN7FXn3PrSB+bn53/zPBaLEYvFvnn99UheRDJXy5a+\nNXGvXtC4MYyo1pAwsxUUFFCQoM1040nyhUCbEq9bF3+spC3AR865XcAuM1sIdAQqTPKlaSQvkh2O\nPBLmzoXTT/fTOANLTwBnudID4Ouvv77a54pnumYpcLSZ5ZpZXWAYMKPUMdOBnmZWy8wOBroCq6sa\njGrkRbJH+/YwcyZcfjnMr3Y9nlSm0iTvnNsLjAWeA1YCU5xzq81stJmNKj5mDTAXWAEsBiY551ZV\nNRjVyItkl5NPhscfh6FD4bXXQkeTmSotoUzom1VSQnnIIfDWW2pRKpJtZs6Eyy6DefP8CF++Ldkl\nlCnxxRfw+efQokXoSEQk1c49F/70J9+9UrtLJVba9JNXjbxIdrvkEvjkE9+L/qWX4LDDQkeUGdIq\nyeumq0h2GzsWPv7Yj+gXLIAmTUJHFH1pM12jm64iAnDddRCLwXnnwZdfho4m+tIqyWskLyJmcMcd\nPh9cdBHs2RM6omhLmySv1a4i8rWcHPj7333CV4vimkmbJK+RvIiUVKcOPPaYb1F85ZVqUVxdSvIi\nkrbq1YMZM/xN2JtuCh1NNKVFdc3OnfCf/6hkSkQO1KQJzJkDp54KzZvDz34WOqJoSYsk/+67cMQR\nfh5ORKS0ww7zLYp79fLtiocMCR1RdKRFktdNVxGpzPe+B7Nn+8VSTZtC376hI4qGtBg7az5eROLR\nsaNvaHbxxfCvf4WOJhrSJslrJC8i8ejdGyZN8v1u1q0LHU36S4vpmi1b/OYBIiLxGDQItm+Hs86C\nl1/2u01J2dIiyRcWQuvWoaMQkSj56U/hgw8gLw8WLlSfm/KkxXTNli1K8iJSdddc42cBBgzwpdhy\noOCbhjgHDRrA++9Dw4YpC0VEMsS+ffDDH/ok//jjUDst5icSK9KbhvznP375shK8iFRHTg48+KDf\neGjMGLU/KC14kt+yBVq1Ch2FiERZ3bowbRosXw7jxoWOJr3EleTNLM/M1pjZWjO7uoLjTjGzPWY2\nON4ANB8vIonQsKFfLPXoo3DXXaGjSR+Vzl6ZWQ4wEegLbAWWmtl059yaMo77AzC3KgGoskZEEqVF\nC5g7F047zT+/8MLQEYUXz0i+C7DOObfJObcHmAIMLOO4XwBPAB9UJQCN5EUkkY46CmbNgiuugIKC\n0NGEF0+SbwVsLvF6S/HHvmFmLYFBzrm7gSrdAVaSF5FEO/FEmDLF7yy1YkXoaMJKVLHRBKDkXH25\niT4/P/+b57FYjC1bYgyOewZfRCQ+p58Of/0rnHMOLFoUrdYpBQUFFCToz5BK6+TNrBuQ75zLK359\nDeCcc7eUOGbD10+B5sAXwCjn3IxS5zqgTv4HP4DJk6FDh5r+U0REDvTnP8Pdd/v2B82ahY6mempS\nJx9Pkq8FvI2/8foe8BpwsXNudTnHPwA845x7sozPHZDkv/tdeOcd3zpURCQZrr4aXnoJXngBDj44\ndDRVl9TFUM65vcBY4DlgJTDFObfazEab2aiyviTeN//8c9i92yd6EZFkGT/e96MfNgyKikJHk1pB\n2xq8/Tacdx6sXZuyEEQkS331lc83ublw771g1RoXhxHZtgaqrBGRVKlbF554ApYtgxtuCB1N6gRt\n5aOWBiKSSg0b+hr6Hj18D/rLLw8dUfIFTfJa7SoiqXbooTBnjt8U/PDD/Q5TmUzTNSKSddq1g+nT\nYeRIWLw4dDTJpSQvIlmpSxffonjQoMwu/lCSF5Gsdc45cOONfgvB998PHU1yBL/xqiQvIiFddpm/\nP3jOOb6hWaZtYBSsTn73bmjUyG/ZlRN86xIRyWbOwahRsHkzPPOM360unUSyTn7rVn9nWwleREIz\n8/1tatf2yT6TthAMlmI1VSMi6aR2bZg6FVauhBLNciMv2Jy8kryIpJv69WHmTL9YqnXrzFgspSQv\nIlJCixbw7LN+sVTLlv6GbJQFna5RSwMRSUft2sHTT/vFUq+/HjqamgmW5NXSQETSWdeu8Le/wYAB\nsGFD5cenK03XiIiUY8AAPyDNy4NXXoHmzUNHVHWarhERqcCYMXDBBb4X/Zdfho6m6oIshtq7F+rV\n8ztD1a2bsrcXEamWfftgxAif5J94AmrVSu37R24x1Pbt0LixEryIRENODvz97/Dpp/CrX0VrsVRc\nSd7M8sxsjZmtNbOry/j8cDN7o/ixyMxOqOh8773nV7uKiETFQQfBk0/C/Pnwpz+FjiZ+ld54NbMc\nYCLQF9gKLDWz6c65NSUO2wD0cs59amZ5wH1At/LOqSQvIlHUpImvoe/RA444AoYODR1R5eKprukC\nrHPObQIwsynAQOCbJO+cK9l2fzFQ4S3V996Dww6rerAiIqEdcYTfQvCMM/xgtVev0BFVLJ7pmlbA\n5hKvt1BxEr8MeLaiE77/vkbyIhJdHTrAI4/AkCGwenXoaCqW0BuvZtYHGAkcMG9fkqZrRCTqzjwT\nbr0V+vf3OS1dxTNdUwi0KfG6dfHHvsXMOgCTgDzn3CflnSw/P58XXvCthjt2jBGLxaoYsohIevjx\nj+Hdd/1m4AsWQIMGiTlvQUEBBQUFCTlXpXXyZlYLeBt/4/U94DXgYufc6hLHtAFeBEaUmp8vfS7n\nnKNnT7j55vSfyxIRqczXG44UFsKMGb5lcaIltU7eObcXGAs8B6wEpjjnVpvZaDMbVXzYdUBT4C4z\nW25mr1V0Ts3Ji0imMIO77vLJ/oor0q+GPuUrXvftczRoANu2Je5PGxGR0D77zM9ODBkC116b2HPX\nZCSf8gZln33mf/MpwYtIJmnY0JdWdu8ObdrAD38YOiIv5UlelTUikqlatoTZs6FPH9+AsU+f0BEF\n6F2j+XgRyWTHHw9TpvjVsCtXho4mQJLXalcRyXSnnw633+5r6LduDRuLpmtERJLghz+ETZt8Df3C\nheHuQwYZySvJi0g2uPZaOOkkP3VTVBQmBs3Ji4gkydc19Hv3wtixYWroNScvIpJEderA44/DkiW+\n102qaU5eRCTJStbQ5+bCsGGpe++UJ3lN14hINmrZEmbOhL59/fNU9e5KeVuDOnUcu3b5PRNFRLLN\nCy/AJZdHpzNpAAAFFUlEQVT4rpXHHRff10RqI+8WLZTgRSR7nXEG3HKLr6Hfti3575fydKupGhHJ\ndpdeCiNGwIAB8OWXyX0vJXkRkQDy8/10zfDhvsQyWVKe5FU+KSLia+jvu8935r3yyuTV0GskLyIS\nSN26MG0azJsHEyYk5z1SXkKpJC8isl+TJr49cY8evg/9BRck9vwpT/KarhER+bY2bfz+sGef7Wvo\nu3dP3Lk1XSMikgY6d4Z//AMGD4b16xN33riSvJnlmdkaM1trZleXc8xfzGydmf3bzE4s71xK8iIi\nZevfH8aN8//dvj0x56w0yZtZDjAROBs4HrjYzI4rdUw/4HvOuXbAaOCe8s536KE1ijdjFBQUhA4h\nbeha7KdrsV+2Xouf/QzOPx8GDoRdu2p+vnhG8l2Adc65Tc65PcAUYGCpYwYCDwE455YAjc2szHR+\n0EE1iDaDZOs3cFl0LfbTtdgvm6/F+PF+j9hLL4V9+2p2rniSfCtgc4nXW4o/VtExhWUcIyIiccjJ\n8fPzW7b4jUdqIuXVNSIiUrnvfAemT695pU2lXSjNrBuQ75zLK359DeCcc7eUOOYeYL5zbmrx6zVA\nb+fctlLnCrAviohI9FW3C2U8I/mlwNFmlgu8BwwDLi51zAzg58DU4l8K/ymd4GsSpIiIVE+lSd45\nt9fMxgLP4efw73fOrTaz0f7TbpJzbraZ9Tez9cAXwMjkhi0iIvFI6aYhIiKSWklZ8ZrIxVNRV9m1\nMLPhZvZG8WORmZ0QIs5UiOf7ovi4U8xsj5kNTmV8qRTnz0jMzJab2VtmNj/VMaZKHD8jjcxsRnGu\neNPMLg0QZtKZ2f1mts3MVlRwTNXzpnMuoQ/8L471QC5QB/g3cFypY/oBs4qfdwUWJzqOdHjEeS26\nAY2Ln+dl87UocdyLwExgcOi4A35fNAZWAq2KXzcPHXfAa/FbYPzX1wHYDtQOHXsSrkVP4ERgRTmf\nr1beTMZIPqGLpyKu0mvhnFvsnPu0+OViMnd9QTzfFwC/AJ4APkhlcCkWz7UYDkxzzhUCOOc+SnGM\nqRLPtXBAw+LnDYHtzrmiFMaYEs65RcAnFRxSrbyZjCSvxVP7xXMtSroMeDapEYVT6bUws5bAIOfc\n3UAmV2LF831xDNDUzOab2VIzG5Gy6FIrnmsxEWhvZluBN4Bfpii2dFOtvKnFUGnCzPrgq5J6ho4l\noAlAyTnZTE70lakNdAZOB+oDr5rZq865BPYnjIyzgeXOudPN7HvA82bWwTn3eejAoiAZSb4QaFPi\ndevij5U+5ohKjskE8VwLzKwDMAnIc85V9OdalMVzLU4GppiZ4ede+5nZHufcjBTFmCrxXIstwEfO\nuV3ALjNbCHTEz19nkniuxUhgPIBz7h0z+z/gOOD1lESYPqqVN5MxXfPN4ikzq4tfPFX6h3QG8CP4\nZkVtmYunMkCl18LM2gDTgBHOuXcCxJgqlV4L59xRxY8j8fPyV2Rggof4fkamAz3NrJaZHYy/0bY6\nxXGmQjzXYhNwBkDxHPQxwIaURpk6Rvl/wVYrbyZ8JO+0eOob8VwL4DqgKXBX8Qh2j3OuS7iokyPO\na/GtL0l5kCkS58/IGjObC6wA9gKTnHOrAoadFHF+X9wIPFiitPA3zrmPA4WcNGY2GYgBzczsXWAc\nUJca5k0thhIRyWAp3/5PRERSR0leRCSDKcmLiGQwJXkRkQymJC8iksGU5EVEMpiSvIhIBlOSFxHJ\nYP8PATF/dZ97zh0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1094dd8d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x_flow_plot, y_flow_plot[1], label=\"$f'$\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Problem with a shock layer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From [1]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equation:\n",
    "$$\n",
    "\\epsilon y'' + x y' = -\\epsilon \\pi^2 \\cos \\pi x - \\pi x \\sin \\pi x\n",
    "$$\n",
    "Boundary conditions:\n",
    "$$\n",
    "y(-1) = -2, y(1) = 0\n",
    "$$\n",
    "\n",
    "I solve it with $\\epsilon=10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "eps = 1e-3\n",
    "\n",
    "def fun_shock(x, y):\n",
    "    t = np.pi * x\n",
    "    return np.vstack((\n",
    "        y[1], -(x * y[1] + eps * np.pi**2 * np.cos(t) + t * np.sin(t)) / eps\n",
    "    ))\n",
    "\n",
    "def bc_shock(ya, yb):\n",
    "    return np.array([ya[0] + 2, yb[0]])\n",
    "\n",
    "x_shock = np.linspace(-1, 1, 10)\n",
    "y_shock = np.zeros((2, x_shock.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Iteration    Max residual    Total nodes    Nodes added  \n",
      "       1          1.36e+00          10             18       \n",
      "       2          2.16e+00          28             44       \n",
      "       3          5.11e-01          72             14       \n",
      "       4          1.20e-02          86             18       \n",
      "       5          3.98e-03          104             3       \n",
      "       6          1.97e-03          107             1       \n",
      "       7          9.28e-04          108             0       \n",
      "Solved in 7 iterations, number of nodes 108, maximum relative residual 9.28e-04.\n"
     ]
    }
   ],
   "source": [
    "res_shock = solve_bvp(fun_shock, bc_shock, x_shock, y_shock, verbose=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_shock_plot = np.linspace(-1, 1, 100)\n",
    "y_shock_plot = res_shock.sol(x_shock_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x10a3e2e80>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8lnP+x/HXp1IUmhgaonNM9n15OCJ0R6kwyggRY+Rn\n7AxjyM+YDsYQP1v2NWHIGinRejO2LBWiVEypkC0hW3U+vz++N87k3Ge7rnNf9/J+Ph7n0X2f++q6\nPp3H3Xnf13c1d0dEREpTs6QLEBGR5CgERERKmEJARKSEKQREREqYQkBEpIQpBERESljkEDCzjcxs\nkpm9ZWZvmtnpWY4bamZzzGy6me0Y9boiIhJdixjOsQI4y92nm9mawGtmNs7dZ/14gJn1Bjq5+2Zm\nthtwM9A5hmuLiEgEke8E3P0jd5+eefw1MBPosMphfYC7M8dMAdqaWfuo1xYRkWhi7RMws3JgR2DK\nKi91ABZUe76IXwaFiIjkWGwhkGkKehg4I3NHICIieS6OPgHMrAUhAO5x98drOGQRsHG15xtlvlfT\nubSYkYhIA7m7NebvxRICwJ3A2+5+bZbXRwGnAA+YWWfgC3dfnO1kWtQuHpWVlVRWViZdRtFoyM9z\n2TKYOxc+/hgWL4YWLaBbN2ivnrCf6P0ZH7NG/f4HYggBM+sCDADeNLNpgAP/C5QB7u63uvuTZra/\nmc0FlgHHRr2uSL769FPo3BlatYLf/Cb84v/mGzjpJOjYEXr0gAEDYKedkq5UJIYQcPfngeb1OO7U\nqNcSyXc//AD9+sEhh8CQIf/92ooV8Oqr8OST8PvfQ9u2MHAgHHUUrLNOMvWKaMZwEUulUkmXUFTq\n+nm6w6mnwtprwz//+cvXW7QIdwgXXQTvvgtXXQUvvQSbbgqnnRa+V0r0/swPlm/t72bm+VaTSH0M\nHQq33w7PPw9rrVX/v/fBB3DddXDbbaHf4O9/h+22a7o6pfiYWaM7hhUCIjFYtAi23x5eew3Kyxt3\njq+/hltugSuugK5dYfBg2HrrWMssSeXl5cyfPz/pMmJRVlbGvHnzfvF9hYBIwoYNg6efhhEjop9r\n2TK44Qb4v/+D3/0O/vEP2GCD6OctVZlfkEmXEYts/5YoIaA+AZEYPP007LdfPOdq0wbOOQdmz4Z1\n1w1NQ5dcAt9+G8/5RapTCIhEtHIlTJgQXwj86Fe/gssvhylTYNo02GYbGDs23muIKAREIpo6NcwF\n2Gijpjl/p07w8MNw001h9NFhh4XOZJE4KAREInr6aejZs+mv07MnzJgBm28OO+wAd94ZhqWKRKEQ\nEIlo3Lj4m4KyWWON0FE8cWIYVnrggWFkkkhjKQREIvjyy9Bev/feub3u9tvDyy9DRUVYfuL++3N7\nfSkeCgGRCCZNgt13h9atc3/t1VYLcwmeegoqK+GPf4Svvsp9HVLYFAIiEeSyKSibnXcOk9SaNw+P\nX3012XqksCgERCLIVadwXdZcE+64I/QX9O4NN96oTuOGMIv+1VgLFy7k0Ucf5YgjjgBg+fLl9OjR\nI6Z/Wd0UAiKNNHdumMC17bZJV/Kzww+HF14Iy08cdVRYikLq5h79q7FmzZpFRUUFH2TG/b744ouU\nN3btkUZQCIg00oQJYW+AKJ8Cm8Jmm8GLL0LLlrDbbmHmseSv7t27M3z4cAYMGADAxIkT2S+HbYwK\nAZFGeuut/N0YpnXrMI/gjDNgzz1D57HkrylTprDnnnsCMGnSJLp3756zaysERBpp9mzYYoukq8jO\nDP70J3j0UTj22LA6qfoJ8lPfvn0ZPXo0Q4cOZcWKFbRr1y5n11YIiDTSO++E2bv5bs89w/pDI0bA\nMcfA998nXZFUN2nSJObOncs555zDkiVLOP3003N6fS0lLdII334L7dqFjtcWkTdpzY1vvoGjj4ZP\nPoGRI8MKpaUg35eSfv3115k6dSqrrbYaLVq0oH///lmPbYqlpGMJATO7AzgQWOzu29fwelfgceC9\nzLcedfd/ZDmXQkDy3owZcOihMHNm0pU0TFUVnHdeaCIaM6Yw7mSiyvcQaIimCIG4PsMMA64D7q7l\nmGfd/aCYrieSqNmzC/MXaLNmMGRIqH3vveGRR6BLl6SrkiTF0ifg7s8BS+o4LM8G0ok0XqGGwI+O\nOw6GD4eDDw5NQ1K6ctkxvLuZTTezMWamnVOloBV6CECY6Tx2LJxySphhLKUpVyHwGtDR3XcErgce\ny9F1RZpEMYQAwC67wHPPwTXXwAUXaAhpKcrJuAZ3/7ra47FmdqOZrePun9d0fGVl5U+PU6kUqVSq\nyWsUaYh8nyPQEL/9bQiC3r3h00/h+uvDYnSSv9LpNOl0OpZzxTZE1MzKgSfcfbsaXmvv7oszjyuA\nB929PMt5NDpI8tqSJVBWBkuX5t+SEVF8+SX06RO2yrz77rDsRDEoLy9n/vz5SZcRi7KyMubNm/eL\n7+fDENH7gBSwLrAYGAy0BNzdbzWzU4CTgOXAt8CZ7j4ly7kUApLXXn4ZTj65OJds/u47OPLIMKfg\n0UeT2SdBGi7xEIiTQkDy3b33wpNPwn33JV1J01ixAgYOhPnz4YknYO21k65I6hIlBLRshEgDFUun\ncDYtWsBdd8HWW0P37vB5jT13UiwUAiINVChrBkXRrFkYNrr33pBKwccfJ12RNJUCWfVEJH8U+53A\nj8zCyqOtW4cgmDgRNtgg6aokbgoBkQZwhzlzwsYtpcAMLroobGqfSsGkSdChQ9JVSZwUAiIN8MEH\nYT/ftm2TriS3LrggDBnt2jUEQceOSVckcVEIiDRAMU0Sa6hzzw1B0K0bTJ6sICgWCgGRBiiV/oBs\nzjwz/KkgKB4KAZEGKPUQAAVBsVEIiDTA7Nlh2GSpqx4EzzwDG22UbD3SeAoBkQYopZFBdTnzzLBT\n2Y9BsOGGSVckjaEQEKkn97CUQllZ0pXkj7/8JSwzsc8+oWlI8wgKj0JApJ4+/jgMD23TJulK8su5\n5/4cBOl0WIVUCodCQKSe5s9XJ2g2558Py5dDjx7hjmDddZOuSOpLaweJ1JOagmo3eDDsvz/stx98\n8UXS1Uh9KQRE6kkhUDszuPRS2Gsv6NULvvoq6YqkPhQCIvWkEKibGVx9Ney8MxxwACxblnRFUheF\ngEg9KQTqxyzsU9ypE/TtG3Yrk/ylEBCpp/ffVwjUV7NmcPvtoYO4Xz/44YekK5JsFAIi9aQ7gYZp\n3hzuuSfsVDZgQBhGKvknlhAwszvMbLGZvVHLMUPNbI6ZTTezHeO4rkiufPllGAK5zjpJV1JYVlsN\nHngg/PwGDgwzjCW/xHUnMAzome1FM+sNdHL3zYATgJtjuq5ITvx4F2CN2sq7tLVqBSNHhp/hKaeE\nmdeSP2IJAXd/DlhSyyF9gLszx04B2pqZ5hVKwVBTUDStW8Po0TB1Kpx9toIgn+SqT6ADsKDa80WZ\n74kUBIVAdGutBU89FfYqHjw46WrkR3m5bERlZeVPj1OpFKlUKrFaREAhEJd27WDcuLBNZZs2Yd0h\nabh0Ok06nY7lXOYx3ZeZWRnwhLtvX8NrNwOT3f2BzPNZQFd3X1zDsR5XTSJxOfzwMOb9iCOSrqQ4\nLFoU9mU480w49dSkqyl8Zoa7N6rHKs7mIMt81WQU8AcAM+sMfFFTAIjkK90JxKtDh9AsdPnlMGxY\n0tWUtliag8zsPiAFrGtm7wODgZaAu/ut7v6kme1vZnOBZcCxcVxXJFcUAvErL4fx48MS1GusAf37\nJ11RaYqtOSguag6SfPPdd9C2LXzzTZgAJfF6882wBPUtt0CfPklXU5iiNAflZcewSD5ZsCA0XygA\nmsZ228GYMdC7N6y+OvTMOuNImoKWjRCpg5qCmt4uu8Bjj8HRR4f9iiV3FAIidVAI5MYee4QlJg49\nFF58MelqSodCQKQOCoHc6dYN7r47DMedOjXpakqDQkCkDgqB3OrVC269NWxV+eabSVdT/NQxLFIH\nhUDu9ekD338fOoknTICtt066ouKlEBCpg0IgGYcdFjaj6dEDJk+GzTdPuqLipBAQqcXKlfDBB7Dx\nxklXUpqOOioEwb77QjodtqyUeCkERGrxwQdhi8RWrZKupHQNHBg29Nlnn3BH8NvfJl1RcVEIiNRi\n3jw1BeWDE04I21Pus0+4IygvT7qi4qEQEKnFnDmw2WZJVyEQdiWrqvr5jkDhHA+FgEgt5sxRh2Q+\nOe20EATduoU7go4dk66o8GmegEgtZs/WnUC+OeMMOP10SKXCyC2JRncCIrVQc1B++vOfoVmzEAST\nJ6uPIAqFgEgWVVUwd65CIF+dfnpY2TWVgkmTNGqosRQCIlksWhT2EVhrraQrkWxOOeXnO4IJE9R/\n0xgKAZEs1BRUGE46Kczj6NYt7FSmJSYaRiEgksXs2fpkWSgGDgxBsO++MHYs7Lhj0hUVDoWASBa6\nEygsAwb8vDPZ449D585JV1QYYhkiama9zGyWmc02s3NreL2rmX1hZlMzX3+L47oiTUlzBArPIYfA\nsGHwu9/BxIlJV1MYIoeAmTUDrgd6AtsAR5jZljUc+qy775z5+kfU64o0Nc0RKEz77w+PPAJHHBHu\nCKR2cdwJVABz3H2+uy8HRgB9ajjOYriWSE6sWBHWDdKqlYVp771D38CJJ8Lw4UlXk9/iCIEOwIJq\nzxdmvreq3c1supmNMTP130tee/99aN8e1lgj6UqksXbZJUwkGzwYrrgi6WryV646hl8DOrr7N2bW\nG3gMyNraWllZ+dPjVCpFKpVq6vpE/ouagorDllvCc8+FLSs/+iiEQbMiWCwnnU6TTqdjOZe5e7QT\nmHUGKt29V+b5IMDdfUgtf+c/wC7u/nkNr3nUmkSiuu46ePttuOmmpCuROCxZEjqLy8vhzjuhZcuk\nK4qXmeHujWpyjyMTXwE2NbMyM2sJ9AdGrVJg+2qPKwjh84sAEMkXGhlUXNq1CxPJvvkmDCFdsiTp\nivJH5BBw95XAqcA44C1ghLvPNLMTzOxPmcP6mdkMM5sGXAMcHvW6Ik1JzUHFZ4014KGHwkSyLl1C\nx7/E0BwUNzUHST7o1AmefBK22CLpSqQpXHstXH45jBwJFRVJVxNdlOYghYDIKn74ISwa99VXxdd2\nLD8bNQqOOy70//Tvn3Q10STdJyBSVN57DzbeWAFQ7A46KMwqHjQI/v73sHR4KVIIiKxCawaVju23\nh5dfDmHQrx98+WXSFeWeQkBkFVo9tLSsv37YlGb99WG33WDWrKQryi2FgMgqZs5UCJSaVq3g5pvh\n7LPDkhOPPZZ0RbmjjmGRVWy+eRhKuMMOSVciSXj5ZTj00NA8dOmlhdE3pI5hkZgsXAiffw7bbZd0\nJZKUigqYOjX0De21F/znP0lX1LQUAiLVTJ4c9qsthvVlpPHWXTcsQ3344aGf4KGHkq6o6eitLlLN\npEmwzz5JVyH5wAzOOgtGj4bzz4ejj4alS5OuKn4KAZEMd4WA/FJFBUybFiYQ7rADxLR4Z95QCIhk\nvPceLF+upSLkl9q0gRtvDF9HHRU2qymWuwKFgEjG5MnhLsC0B55ksf/+8NZb4a5x223hiSeSrig6\nhYBIhpqCpD7atoVbboF77gl9Bn37FvYIIoWACD/3B3TrlnQlUihSKZgxI4we2nVXqKwM+xUUGoWA\nCGGW8BprwCabJF2JFJJWreC880LH8cyZoT/pjjtgxYqkK6s/hYAIagqSaDbeGB54IMwnuOeesDDd\no48WxsqkCgERfu4UFomic+fwXrrySrjkkjCk9L778vvOQGsHSclbtix8kpsxAzbcMOlqpFi4w9NP\nhzD46CM4/XQ45hhYe+34r6W1g0QiuPBC6N1bASDxMoNeveDf/4Zhw+C556C8HE4+GV5/PenqfhZL\nCJhZLzObZWazzezcLMcMNbM5ZjbdzHaM47oiUb3+Otx1F1x9ddKVSDHbc8/QZzBjRti34KCDwiKF\nQ4bA++8nW1vk5iAzawbMBvYFPgBeAfq7+6xqx/QGTnX3A8xsN+Bad++c5XxqDpKcWLkS9tgDjj8e\n/ud/kq5GSklVFTz/PNx7Lzz8MJSVwQEHhK9dd4XmzRt2vkQ3mjezzsBgd++deT4IcHcfUu2Ym4HJ\n7v5A5vlMIOXui2s4n0JAcuKGG8Kns3Raq4ZKclasgBdfDAvVjRkDCxaEuQdduoSO5m22gQ4dap/J\nHiUEWjS28Go6AAuqPV8IVNRxzKLM934RAiK5MG9emNzzzDMKAElWixZh34K99grNQ59+Ci+8EPoQ\nLrsszD/49lvYaqswgGHDDcNXu3Zhbsvqq0e8fjz/jHhVVlb+9DiVSpFKpRKrRYrLypVhEbCLLgoh\nsPXWSVck8t9+/evQZ3DQQT9/77PPwt7HCxfCBx/AlClp5sxJs3x59OGncTUHVbp7r8zz+jQHzQK6\nqjlIcsU9fLo6/fQwRO/GG8MnK5FikPQQ0VeATc2szMxaAv2BUascMwr4A/wUGl/UFAAicfvsM7j2\n2jBp5+ij4YwzwuxgBYBIELk5yN1XmtmpwDhCqNzh7jPN7ITwst/q7k+a2f5mNhdYBhwb9boi2VRV\nhVmbt90GTz0VRlxcey107ar2f5FVacawFI2lS+H22+Gmm6B16zD086ijQgeaSDFLenSQSKLefz98\n0r/rLujZM4y93m03bQ4jUh+6OZaCtWBB2OZvp53C82nTwmJdnTsrAETqSyEgBeejj+DUU0Nn769+\nBe+8E1Zt7Ngx6cpECo9CQArGt9/CP/8Z9nZt2TKMm77ssjCuWkQaR30Ckvfcw/oqf/0r7LILTJkC\nnTolXZVIcVAISF577z045ZQwU3L48DDMU0Tio+YgyUvLl8Oll0JFRdjQe+pUBYBIU9CdgOSdGTPC\nDkzrrQevvKLN30Waku4EJG+sWBE+/XfrBiedBGPHKgBEmpruBCQvzJsHAwaEpXFffTVssiEiTU93\nApK4Bx4Ibf8HHwzjxikARHJJdwKSmG++gdNOCxtxjx0bhn+KSG7pTkASMWcO7L47fPddGPmjABBJ\nhkJAcm7kyLB/6oknhsXe1lwz6YpESpeagyRnVq6ECy4Ii7yNGQO77pp0RSKiEJCcWLo0jP75+usw\n9n+99ZKuSERAzUGSA7Nnh+Wdy8th/HgFgEg+UQhIk0qnYa+94Kyz4PrrYbXVkq5IRKpTc5A0mWHD\nYNAguP9+2GefpKsRkZpECgEzawc8AJQB84DD3H1pDcfNA5YCVcByd6+Icl3Jb1VVoQN4xAh45hnY\ncsukKxKRbKI2Bw0CJrj7FsAk4Lwsx1UBKXffSQFQ3H74ISz+NmkSvPSSAkAk30UNgT7A8Mzj4UDf\nLMdZDNeSPPfll3DAAeHPiRPVASxSCKL+Yl7f3RcDuPtHwPpZjnNgvJm9YmbHR7ym5KEPPwzr/W+6\nKTzyCLRunXRFIlIfdfYJmNl4oH31bxF+qf+thsM9y2m6uPuHZrYeIQxmuvtz2a5ZWVn50+NUKkUq\nlaqrTEnQu+/CfvvBscfC+eeDWdIViRS3dDpNOp2O5Vzmnu33dj3+stlMQlv/YjP7DTDZ3beq4+8M\nBr5y96uyvO5RapLceuMN6N0b/v53OOGEpKsRKU1mhrs36uNX1OagUcAfM4+PAR5f9QAza21ma2Ye\ntwH2A2ZEvK7kgeefhx494OqrFQAihSrqncA6wIPAxsB8whDRL8xsA+A2dz/QzDYBRhKailoA/3L3\ny2o5p+4ECsD48XDkkWEBuJ49k65GpLRFuROIFAJNQSGQ/x5/HI4/Hh59FPbcM+lqRCTJ5iApMfff\nH5p+xo5VAIgUA4WA1NuwYXD22TBhgjaBESkWWjtI6uXWW+Hii8NM4C22SLoaEYmLQkDqdMMNcPnl\nMHlymAwmIsVDISC1Gjo0DAFNp2GTTZKuRkTiphCQrK69Nnyl01BWlnQ1ItIU1DEsNfoxACZPVgCI\nFDOFgPyCAkCkdKg5SP7LddfBNdeoCUikVCgE5Cc33QRXXqkAECklCgEB4Lbb4NJLQxNQeXnS1YhI\nrigEhOHD4cILQwB06pR0NSKSSwqBEnf//XDeeWEm8GabJV2NiOSaQqCEPfIInHlmWAtIG8KLlCaF\nQIkaPRpOPhmeegq23TbpakQkKQqBEjR+PAwcGIJgp52SrkZEkqQQKDHPPht2BBs5Eioqkq5GRJKm\nGcMlZMoUOOSQ0BmsDWFEBBQCJWP6dDjoILjrLujePelqRCRfRAoBM+tnZjPMbKWZ7VzLcb3MbJaZ\nzTazc6NcUxpu5kzo3TvsC3DAAUlXIyL5JOqdwJvAwcAz2Q4ws2bA9UBPYBvgCDPTgMQcefdd2G+/\nsClMv35JVyMi+SZSx7C7vwNgZrXtcl8BzHH3+ZljRwB9gFlRri11W7AgNP387W9w9NFJVyMi+SgX\nfQIdgAXVni/MfE+a0Ecfwb77wmmnwQknJF2NiOSrOu8EzGw80L76twAHznf3J5qiqMrKyp8ep1Ip\nUqlUU1ymaH32GfToET79n3VW0tWISNzS6TTpdDqWc5m7Rz+J2WTgL+4+tYbXOgOV7t4r83wQ4O4+\nJMu5PI6aStXSpeEOoHv3sCporQ11IlIUzAx3b9T/9jibg7IV8AqwqZmVmVlLoD8wKsbrSsayZWH0\nT+fOCgARqZ+oQ0T7mtkCoDMw2szGZr6/gZmNBnD3lcCpwDjgLWCEu8+MVras6rvvoE8f2HxzGDpU\nASAi9RNLc1Cc1BzUcD/8AL//Pay1Ftx7LzRvnnRFIpJLUZqDFAIFbsUK6N8//PnQQ7DaaklXJCK5\nFiUEtIBcAVu5Eo49Fr76CkaNUgCISMMpBApUVRWceCIsXAhjxkCrVklXJCKFSCFQgNzhjDPg7bfh\n6aehdeukKxKRQqUQKDDucM458NJLYVvINddMuiIRKWQKgQLiHtYBGj8+bAzftm3SFYlIoVMIFJCL\nL4bHH4fJk2GddZKuRkSKgUKgQFx2WdgRLJ2G9dZLuhoRKRYKgQJwxRVwxx3wzDPQvn3dx4uI1JdC\nIM9ddRXccku4A9hww6SrEZFioz2G89g118CNN4Y+gI02SroaESlGuhPIU9deGxaCS6dh442TrkZE\nipVCIA9ddVXYFD6dho4dk65GRIqZQiDPXHHFz30AugMQkaamEMgjl10Gd94ZAkB9ACKSCwqBPOAO\nF14II0aETuAOHZKuSERKhUIgYe4waBCMHat5ACKSewqBBFVVwZ//DM8/H+4A1l036YpEpNRE3WO4\nn5nNMLOVZrZzLcfNM7PXzWyamb0c5ZrFYsUKOO44ePVVmDhRASAiyYh6J/AmcDBwSx3HVQEpd18S\n8XpF4fvvYcAAWLoUxo3TctAikpxIIeDu7wCYWV17WxqanQzAsmVhU/g2bWD0aO0IJiLJytUvZgfG\nm9krZnZ8jq6Zdz77DLp3hw02gAcfVACISPLqvBMws/FA9TErRvilfr67P1HP63Rx9w/NbD1CGMx0\n9+caXm7hWrAAevaEAw+EIUOgznsnEZEcqDME3L1H1Iu4+4eZPz8xs5FABZA1BCorK396nEqlSKVS\nUUtI1NtvQ+/eYV/gs85KuhoRKXTpdJp0Oh3Luczdo5/EbDJwtru/VsNrrYFm7v61mbUBxgEXuvu4\nLOfyOGrKF88+C4ceCldeCUcdlXQ1IlKMzAx3b1T7QtQhon3NbAHQGRhtZmMz39/AzEZnDmsPPGdm\n04CXgCeyBUCxuf9+6NcP/vUvBYCI5KdY7gTiVAx3Au5hHaCbbw4jgLbbLumKRKSYRbkT0IzhmH3/\nPZx0EkyfDi++qN3ARCS/aex+jD75JAwB/eIL+Pe/FQAikv8UAjGZMQMqKqBrV3j44TAZTEQk36k5\nKAYPPginnBL2BB4wIOlqRETqTyEQwYoVcN554ZP/uHGw005JVyQi0jAKgUZavBiOPBKaNw8rgWoV\nUBEpROoTaIRJk2DnnWGPPcJmMAoAESlUuhNogJUr4eKL4dZbYfhw6BF5QQ0RkWQpBOrpP/+BP/wB\nWraE114LK4GKiBQ6NQfVwR3uuisM/+zbF8aPVwCISPHQnUAtPvwwDP2cOzdsAbn99klXJCISL90J\n1MAd7rwTdtgBttoKXn5ZASAixUl3AquYPTt8+v/88zD2f8cdk65IRKTp6E4g4+uvw8SvPfaAXr1g\nyhQFgIgUv5IPgaqqsN7/VlvBwoXwxhvwl79AC90jiUgJKNlfde6huWfQoDDs8777YK+9kq5KRCS3\nSi4E3CGdDpO+Fi2CSy+Fgw/Wxu8iUppKJgSqquDJJ+GSS0Kn76BBcPTRavYRkdJW9L8Cly4NSzzc\ncAOssQb87//CIYeEhd9EREpd1I3mLzezmWY23cweMbO1sxzXy8xmmdlsMzs3yjXro6oKnn0Wjj8e\nysvhhRfg9tth2jQ47DAFgIjIj6KODhoHbOPuOwJzgPNWPcDMmgHXAz2BbYAjzGzLiNf9hZUr4aWX\nwjDPTTYJY/033RTeegtGjAidvqXW7p9Op5Muoajo5xkv/TzzQ6QQcPcJ7l6VefoSsFENh1UAc9x9\nvrsvB0YAfaJcN1wb3n0X7r4bjjkmrOdz/PHhtSeegDffhHPPLe19fvWfLF76ecZLP8/8EGefwEDC\nL/hVdQAWVHu+kBAM9eIOS5bA++/D22+Hrxkzwqf+5s2hSxdIpeCii6CsLNo/QESk1NQZAmY2Hmhf\n/VuAA+e7+xOZY84Hlrv7fXEU1aULfPddGMXz4Yew+uqw0Uaw9dbha8CAsJ9vWVnpNfGIiMTJ3D3a\nCcz+CBwP7OPu39fwemeg0t17ZZ4PAtzdh2Q5X7SCRERKkLs36iNxpOYgM+sF/BXYu6YAyHgF2NTM\nyoAPgf7AEdnO2dh/iIiINFzU0UHXAWsC481sqpndCGBmG5jZaAB3XwmcShhJ9BYwwt1nRryuiIjE\nIHJzkIiIFK5EVxE1s35mNsPMVprZzrUcl9PJZoXKzNqZ2Tgze8fMnjaztlmOm2dmr5vZNDN7Odd1\n5rv6vN/MbKiZzclMlNSi41nU9bM0s65m9kWmJWGqmf0tiToLhZndYWaLzeyNWo5p0Hsz6aWk3wQO\nBp7JdkBWzwZVAAACaElEQVSuJpsViUHABHffAphEDZP3MqqAlLvv5O71Hq5bCurzfjOz3kAnd98M\nOAG4OeeFFoAG/N991t13znz9I6dFFp5hhJ9njRrz3kw0BNz9HXefQxh2mk2TTDYrUn2A4ZnHw4G+\nWY4zkv8AkK/q837rA9wN4O5TgLZm1h5ZVX3/72owSD25+3PAkloOafB7sxB+EdQ02axDQrXku/Xd\nfTGAu38ErJ/lOCd05r9iZsfnrLrCUJ/326rHLKrhGKn//93dM00XY8xs69yUVrQa/N5s8lVE6zPZ\nTOqvlp9nTW2p2Xr9u7j7h2a2HiEMZmY+YYjk2mtAR3f/JtOU8RiwecI1lZQmDwF37xHxFIuAjtWe\nb5T5Xkmq7eeZ6TBq7+6Lzew3wMdZzvFh5s9PzGwk4bZdIRDU5/22CNi4jmOkHj9Ld/+62uOxZnaj\nma3j7p/nqMZi0+D3Zj41B2VrF/xpspmZtSRMNhuVu7IKyijgj5nHxwCPr3qAmbU2szUzj9sA+wEz\nclVgAajP+20U8Af4aUb8Fz82w8l/qfNnWb292swqCMPWFQC1M7L/vmzwezPRTWXMrC9hwtmvgdFm\nNt3de5vZBsBt7n6gu680sx8nmzUD7tBks6yGAA+a2UBgPnAYhMl7ZH6ehKakkZnlOVoA/3L3cUkV\nnG+yvd/M7ITwst/q7k+a2f5mNhdYBhybZM35qj4/S6CfmZ0ELAe+BQ5PruL8Z2b3ASlgXTN7HxgM\ntCTCe1OTxURESlg+NQeJiEiOKQREREqYQkBEpIQpBERESphCQESkhCkERERKmEJARKSEKQRERErY\n/wM5QGqQNUxAtwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a3e26a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x_shock_plot, y_shock_plot[0], label='$y$')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Problem with a singular term"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example 6 from [2]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equation:\n",
    "$$\n",
    "y'' + \\frac{2}{x} y' = \\phi^2 y \\exp \\left(\\frac{\\gamma \\beta (1 - y)}{1 + \\beta (1 - y)} \\right)\n",
    "$$\n",
    "Boundary conditions:\n",
    "$$\n",
    "y'(0) = 0, y(1) = 1\n",
    "$$\n",
    "Parameters:\n",
    "$$\n",
    "\\phi = 0.6, \\gamma = 40, \\beta = 0.2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 3 different solutions, we can obtain all of them using different starting guesses. Note that the singular term is passed separately to the solver as a matrix `S`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "phi = 0.6\n",
    "gamma = 40\n",
    "beta = 0.2\n",
    "\n",
    "def fun_sing(x, y):\n",
    "    return np.vstack((\n",
    "        y[1], \n",
    "        phi**2 * y[0] * np.exp(gamma * beta * (1 - y[0]) / ( 1 + beta * (1 - y[0]) ) )\n",
    "    ))\n",
    "\n",
    "def bc_sing(ya, yb):\n",
    "    return np.array([ya[1], yb[0] - 1])\n",
    "\n",
    "S = np.array([[0.0, 0], \n",
    "              [0, -2]])\n",
    "\n",
    "x_sing = np.linspace(0, 1, 10)\n",
    "\n",
    "y_sing_1 = np.empty((2, 10))\n",
    "y_sing_1[0] = 1\n",
    "y_sing_1[1] = 0\n",
    "\n",
    "y_sing_2 = np.empty((2, 10))\n",
    "y_sing_2[0] = 0.5\n",
    "y_sing_2[1] = 0\n",
    "\n",
    "y_sing_3 = np.empty((2, 10))\n",
    "y_sing_3[0] = 0\n",
    "y_sing_3[1] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Iteration    Max residual    Total nodes    Nodes added  \n",
      "       1          8.74e-06          10              0       \n",
      "Solved in 1 iterations, number of nodes 10, maximum relative residual 8.74e-06.\n"
     ]
    }
   ],
   "source": [
    "res_sing_1 = solve_bvp(fun_sing, bc_sing, x_sing, y_sing_1, S=S, verbose=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Iteration    Max residual    Total nodes    Nodes added  \n",
      "       1          1.95e+00          10             14       \n",
      "       2          8.11e-05          24              0       \n",
      "Solved in 2 iterations, number of nodes 24, maximum relative residual 8.11e-05.\n"
     ]
    }
   ],
   "source": [
    "res_sing_2 = solve_bvp(fun_sing, bc_sing, x_sing, y_sing_2, S=S, verbose=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Iteration    Max residual    Total nodes    Nodes added  \n",
      "       1          2.30e-01          10              9       \n",
      "       2          1.72e-03          19              5       \n",
      "       3          9.77e-04          24              0       \n",
      "Solved in 3 iterations, number of nodes 24, maximum relative residual 9.77e-04.\n"
     ]
    }
   ],
   "source": [
    "res_sing_3 = solve_bvp(fun_sing, bc_sing, x_sing, y_sing_3, S=S, verbose=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_sing_plot = np.linspace(0, 1, 100)\n",
    "y_sing_plot_1 = res_sing_1.sol(x_sing_plot)\n",
    "y_sing_plot_2 = res_sing_2.sol(x_sing_plot)\n",
    "y_sing_plot_3 = res_sing_3.sol(x_sing_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x106eca160>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEACAYAAABWLgY0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FNX+x/H3SQihVyE0ERHpTREE4QcBVGJFLJduFyzc\niyISwItERLgXFZUriIAoIEgVBBGlLr2EGgIJVZASQkuAAAkp5/fHZN1NTMiS7O5s+b6e5zwzszvZ\nfBnCh8mZM2eU1hohhBC+KcDsAoQQQriOhLwQQvgwCXkhhPBhEvJCCOHDJOSFEMKHScgLIYQPyzPk\nlVLfKqXilVJRN9lnnFLqkFJqt1KqqXNLFEIIkV+OnMl/B3TK7U2l1CPAXVrru4G+wEQn1SaEEKKA\n8gx5rfUGIOEmu3QGpmfuuxUorZQKcU55QgghCsIZffJVgRN226cyXxNCCGEyufAqhBA+rJATPuMU\ncLvddrXM1/5GKSUT5QghRD5orVV+vs7RkFeZLSeLgbeAOUqplkCi1jo+tw+SCdEMERERREREmF2G\nR5BjYSPHwsZTj0VyMsTFGe306azr1hYXB1euQOXKRqtSxWiVK0PVqrb12yqlsDNhFYtiF7Lk4BIq\nFq9I5zqd6Vy3M80KVUfdcw9Mn4568MF815tnyCulZgGhQHml1J/AcKAwoLXWk7TWvyqlHlVKHQau\nAi/luxohhDBJUpItsHNq1vBOSoJKlbIGd+XK0LatbbtKFShfHgJy6BC/knKFZYeXMTXmJ35b+huN\nQhrxVJ2nCG8TTq1ytYydtIYnn4RevaBjxwL9ufIMea11Dwf26VegKoQQwgW0hosX4cyZv4d29tdS\nU7OeeVvX69XL+nq5cjmH981cvH6RxQcW81PMT1iOWWhTvQ1d6nbhy7AvCSmRw2DECROMohYsKPAx\ncEafvMiH0NBQs0vwGHIsbORY2NzsWKSkQHy8LaytzT68ra8VK2aceVuD2tplct99WV8rXRpUvnq9\ncxafFM+i2EXMj5nPtlPbeLDmg3Rt0JUZXWZQukjp3L9w/36IiICNG6Fw4QLXodzZR66U0tInL4TI\nSXo6nD+fNbTj47NuW1tSEoSEGM0a0tYgr1Qpayta1H1/htNXTvNTzE/M3z+f3Wd28+jdj/JMvWcI\nqxVG8cLF8/6AlBRo2RLefBNee+2vl5VS+b7wKiEvhHCZjAy4cMEW2NbQtl9a1y9cgLJljWAOCbGF\ntDXI7bfz02XiKqcun2JBzALm7Z/HvrP7eLz24zxb/1kevuthihQqcmsfFh4OBw7AwoVZfq2QkBdC\nuE1qKpw7Z4Tz2bO2oM7ezp41zsxLl7addVvD2/4s3PpahQpQyEs6kOOuxLEgZgFz980l+mw0T9R5\ngufqP8dDNR8iuFBw/j7UYoGePWH3buNg2JGQF0Lkm9Zw6ZIRyrk1+0C/csUYOWIN6ooVs4a4fatQ\nAYKCzP4TOse5q+dYELOAOfvmsPvMbh6v/ThdG3QtWLBbJSZCkybwzTcQFva3tyXkhRB/0doI4nPn\njHb2bNZ167b960WKGGFdoYItuO3XrYHtaV0lrpaYnMjCmIXM3jebrSe38sjdj9C1QVfCaoXdelfM\nzfTqZfRV/e9/Ob5dkJD3kl+OhPBfaWnGMMDz54127pxtmX3d2oKCjFC2hrV1vWpVaNrU9rp1WcSJ\neeXtrt64ypKDS/gx+kcsxyx0vLMjr9zzCgu7LqRYUDHnf8N582D7dti50/mfjZzJC+FW6enGb+bn\nzxsXGq3BbV23f80a4JcuQZkycNttRrOGdoUKWbft1905osQXpKansuLoCmbuncnSg0tpWa0l3Rt2\n56m6T918uGNBxcUZ/+v+8gs0b57rbtJdI4SbaQ3XrhmhfOGCcaZtv8y+bg3vS5eMC5HlyxvNGtz2\n69Zta3CXLQuBgWb/iX2P1ppNJzYxc+9M5u2fR+3ytenRsAfPNXiOisUruqMAeOwxaNHCGBd/E9Jd\nI0Q+Wc+sExJs7eJF29Lasm9fvGj0S5cvb4SwNbTLlTOWVapAo0a2162tbFnvGUHiq2LPx/JD1A/M\n2juLIoWK0KtxL7a9uo07y97p3kK++cb4de399136beRMXng1reH6deMMOTHRFth5La3tyhUoVcoI\n33Llsi6t4V2unO0163bZstIl4k3OXj3Lj3t/5Ie9P3Dq8im6NexG78a9aVqpKcqZt7k66sgR46an\n9euhbt08d5fuGuGVtDZm9Lt82QjpS5eyrttvJybaXrOGuXU9MNDosy5d2gjf7Ov2y+ytVCnpCvFV\nyWnJLD6wmOl7prPhzw08WedJejXuRcc7OxIYYOJfeno6hIbCM8/A22879CUS8sJtrH3RSUnGWfCV\nK8b65cu27ezt8mXb+9Z1a4AHBhqBXKqUbWldt2/W4LZfty5lZIiwsvazT9szjQUxC2hWuRm9G/em\nS70ulChcwuzyDGPHws8/w5o1Do9FlZAXf5OaaoTxtWtw9arR7NftW1LS37dzaleuGO8HB0OJElCy\npG2ZUytVyra0rtsHecmSxmcJUVDHE48zfc90pu2ZRlBgEC80eYFejXtRrVQ1s0vLKibGmJN461ao\nWdPhL5OQ92BpacacQ9aWnJx1Pbd2/bptaW3W7WvXcl5ev24L84wMKF7caEWL2tZzayVKZF23bluD\n3BrmxYvLhUPhGa6lXuOnmJ/4bvd37D6zm24NuvFC0xdoXqW5Of3seUlLgwcegFdegb59b+lLvWp0\nzd69xtKa9Vr/fd3aMjJsS2uz305PN5r9em4tLc22zN5SU/++zN5u3Pj7MqeWkmJbpqQYf67gYFsr\nUsS2XrSosW19zX7bfr18eWM7p2YN8WLFsm4HBTl32lQhPIHWmm2ntjF111Tm7Z/H/dXup2+zvjxZ\n50nn3oHqCmPGGH2Mffq49du6/Uy+YUNtt21bZl8PCLCtBwbmvB0QYKxbt63rObVChYxmXQ8Ksm0H\nBRnNft2+FS6cdRkUZIRy4cK21+y3revBwXLWK4QznL16lhl7ZjB191RupN/g5aYv83yT56laqqrZ\npTkmOhrat4cdO6B69Vv+cumuEUL4nPSMdJYfWc6UXVNYdXQVT9V9ilfueYU21dt4ZndMblJToVUr\no4vGbo74W+FV3TVCCHEzxxOPM3XXVKbunkrlEpV55Z5X+K7zd5QKLmV2afkzZozR5/rqq6Z8ezmT\nF0KYLjU9lSUHlzB552S2ndpGj4Y9ePXeV2lSqYnZpRVMAbtprKS7Rgjhlf5I+IPJOyfz3e7vqFWu\nFn3u7cOz9Z+laJAP3E6clmZ007z2WoEvtkp3jRDCa6RlpLHkwBIm7pjIjtM76N24N6ufX029CvXM\nLs25xo41bgzJZz+8s8iZvBDCLU5ePsnkHZOZsmsKd5a5k77N+vJcg+c8f+hjfhw4AK1bQ2Qk3Fnw\nic/kTF4I4ZEydAYrj67k6+1fs/bYWro37M5vPX+jUUgjs0tznYwM44anDz5wSsAXlIS8EMLpEq4n\n8N3u7/h6+9cUDyrOm83fZEaXGZ4zf4wrjR9vLPv1M7eOTNJdI4Rwml1xuxgfOZ4FMQt47O7HeLP5\nm7Sq1sq7xrUXxPHj0KwZbNjg0BTCjpLuGiGEaW6k32DB/gV8FfkVJy6d4PX7XudAvwPuebqSJ9Ha\nuOFpwACnBnxBScgLIfIl7koc3+z4hkk7JlH3troMbDWQJ+o8QaEAP42VmTONZ7a+957ZlWThp38b\nQoj82nZqG+O2jmPpoaV0a9CNFb1X0KBiA7PLMte5czBwoPFA7qAgs6vJQvrkhRB5Sk1PZf7++Xy5\n9Uvir8bTr3k/Xr7nZcoWLWt2aZ6hZ0+oXBk+/dQlHy998kIIl7hw7QKTdkxifOR47i5/N4PbDOaJ\n2k+Y+/g8T7NsGWzZYptH3cNIyAsh/ibmXAxfbPmCufvn0qVuF5b2WOr988i4wtWr8MYbMHmy8VAH\nDyQhL4QAjAdyrPpjFWM3j2Vn3E7euO8NYt+KJaREiNmlea7hw43H+T30kNmV5Er65IXwcylpKcyO\nns3YLWNJz0hnQKsB9GjUwzenG3CmnTvh0UeNbpoKFVz6raRPXghxyxKuJzBx+0S+ivyKBhUaMObB\nMTx818P+c+NSQaSlGROPjRnj8oAvKAl5IfzMscRjfLHlC6bvmc4TdZ5gWc9lNA5pbHZZ3mXcOChb\nFnr3NruSPAU4spNSKkwpFauUOqiUCs/h/VJKqcVKqd1Kqb1KqRedXqkQokB2xu2k+4LuNJvUjODA\nYKLeiGLaU9Mk4G/Vn3/CqFEwcaLt4dQeLM8+eaVUAHAQ6AicBiKBblrrWLt9hgCltNZDlFK3AQeA\nEK11WrbPkj55IdxIa83KoysZs2kMMedieLvl2/Rp1sd7H6VnNq2hc2do0QL+/W+3fVtX98m3AA5p\nrY9nfrPZQGcg1m4fDZTMXC8JXMge8EII90nPSGf+/vmM2TSG66nXGdR6ED0a9aBwYGGzS/NuixbB\noUMwb57ZlTjMkZCvCpyw2z6JEfz2vgIWK6VOAyWArs4pTwhxK5LTkpm2exqfbPqEkBIhDG83nMdr\nP06AcqhnVtzMlSvwr3/BDz9AcLDZ1TjMWRdeOwG7tNYdlFJ3ASuUUo211knZd4yIiPhrPTQ0lNDQ\nUCeVIIT/upJyhYnbJ/L5ls+5p/I9fP/U97Sp3sbssnzLBx/Agw9Cu3Yu/1YWiwWLxeKUz3KkT74l\nEKG1DsvcHgxorfV/7fb5BRittd6Yub0KCNdab8/2WdInL4QTXbh2gXFbxzFh+wQerPkgg1sPljtT\nXWHXLggLg3374Lbb3P7tXd0nHwnUUkrdAcQB3YDu2fY5DjwIbFRKhQC1gaP5KUgIkbe4K3F8tvkz\nvtv9Hc/Ue4bNr2ymVrlaZpflmzIyjKkLRo0yJeALKs+Q11qnK6X6Acsxhlx+q7WOUUr1Nd7Wk4CR\nwPdKqajMLxuktb7osqqF8FPHEo8xZuMY5uybw/ONnyfq9Siqlqpqdlm+bcoUCAyEl14yu5J8kWkN\nhPAChy4cYvSG0fx84Gf63NuHd1q9439PXjLD2bPQsCGsXAmNzbufQKY1EMJHxZyLYeT6kSw/spx+\nzftx+J+HZQ53dxo0yLir1cSALygJeSE8UPTZaEauG8nqP1bzdsu3+fqxr+UGJndbvx5WrYL9+82u\npEAk5IXwIFHxUYxYO4L1f67n3VbvMuXJKZQoXMLssvxPaiq8+SZ89hmULJn3/h5MQl4ID2AN9w1/\nbmDgAwOZ9tQ0ihcubnZZ/mv8eKhUCZ57zuxKCkwuvAphouiz0URYItjw5wbee+A9Xr/vdQl3s8XF\nQaNGsGED1K1rdjVAwS68SsgLYYL95/bz4doPWXtsLQMfGMgb970h4e4pevaE6tVh9GizK/mLjK4R\nwkscOH+AEetGsPLoSga0HMDUJ6dKuHsSi8U4g/fyi632JOSFcIOjCUcZsXYESw8t5e3732biYxMp\nGezdF/R8Tmoq9OsHY8dCcd/5j1dCXggXOnHpBCPXjWR+zHz6Ne/HoX8eokyRMmaXJXIyfjxUqQJP\nP212JU4lIS+EC8QnxTNq/ShmRM2gT7M+HOx3kPLFyptdlshNfDx8/LExNt4LnvZ0KyTkhXCihOsJ\nfLLpE77Z8Q29GvVi/1v7qVSiktllibwMHmzMTeMho2mcSUJeCCdIupHEuK3j+HzL53Su05ldfXdR\nvXR1s8sSjti8GVasgJgYsytxCQl5IQogJS2FSTsmMWrDKNrd0Y6NL2+kdvnaZpclHJWeblxsHTPG\n6+9szY2EvBD5kJ6Rzqy9s/jA8gH1K9Tn1x6/ck/le8wuS9yqqVOhWDHonv0RGb5DboYS4hZorfnl\n4C8MXT2UUsGlGN1xNG3vaGt2WSI/EhKgXj347Tdo2tTsam5K7ngVwg02/LmBwSsHk5icyKiOo3ii\n9hMoHxuJ4Vf694eUFJg40exK8iR3vArhQvvO7mPIqiHsid/DiNAR9Grci8CAQLPLEgWxbx/8+KNP\n3dmaGwl5IXJx4tIJhluG88vBXxjSZghzn5tLkUJFzC5LFJTWxln8sGFe+czWWyUhL0Q2icmJjF4/\nmim7pvB6s9c59M9DlC5S2uyyhLMsXAhnzhgP5/YDEvJCZEpJS2FC5ARGbxhN5zqd5SHZvig5GQYO\nhMmToZB/xJ9//CmFuIkMncGc6DkMXT2UhhUbYnnRQv0K9c0uS7jC558bI2k6djS7EreR0TXCr1mO\nWXhvxXsoFJ889AntarQzuyThKqdPGw/k3rYNatY0u5pbIkMohbhF+8/tJ3xlOPvO7mNUx1H8o8E/\nCFABZpclXOnFF6FyZY96GIijZAilEA6KT4pnuGU4P8X8xOA2g5n/3HyCCwWbXZZwtW3bYPlyOHDA\n7ErcTk5dhF+4lnqNketGUn9CfYoFFSO2XywDWg2QgPcH1iGTH3/ss/PT3IycyQuflqEz+CHqB95f\n/T4tq7Vk26vbuKvcXWaXJdzpxx+Npz698ILZlZhCQl74LMsxC+8uf5eggCDmPDuHB25/wOyShLtd\nu2bMFT9rFgT4Z8eFhLzwOYcuHGLQykHsitvFfx78D10bdJU5ZvzVp59Cq1bQpo3ZlZhGRtcIn5Fw\nPYERa0cwI2oG7z3wHv1b9pdpCPzZyZPQpAns2AE1aphdTYEUZHSNf/7+InxKanoq47aOo85XdUhO\nS2b/W/sJbxMuAe/vhg6Fvn29PuALSrprhNfSWrP00FIGLh9I9dLVWfX8KhqFNDK7LOEJIiNh5Uq/\nHDKZnYS88ErRZ6MZ8PsATlw+wdhOY3mk1iPS7y4MWsM778BHH/nlkMnspLtGeJVzV8/x5tI36TCt\nA0/WeZKo16N49O5HJeCFzYIFkJRk3OEqJOSFd7iRfoOxm8dSf0J9CgcWJrZfLP1a9CMoMMjs0oQn\nSUmBQYPgs88gUB7sAtJdIzyc1polB5fw7vJ3qVO+DutfWk/d2+qaXZbwVOPGQcOGfjXLZF4cGkKp\nlAoDvsA48/9Wa/3fHPYJBT4HgoBzWuv2OewjQyiFw6LPRvPO7+9w6vIpPu/0OZ1qdTK7JOHJzp0z\nHsy9cSPUqWN2NU7l0lkolVIBwEGgI3AaiAS6aa1j7fYpDWwCHtZan1JK3aa1Pp/DZ0nIizydv3ae\nD9Z8wPz98/mg3Qf0bdZXumVE3t56y+iiGTfO7EqcztWzULYADmmtj2d+s9lAZyDWbp8ewAKt9SmA\nnAJeiLykpqcyIXICI9ePpFuDbsS8FUP5YuXNLkt4g5gYmDsXYmPz3tfPOBLyVYETdtsnMYLfXm0g\nSCm1BigBjNNaz3BOicIfLDu0jAHLB1C9dHUsL1hoULGB2SUJbzJoEISHQ3k5KcjOWRdeCwH3Ah2A\n4sBmpdRmrfVhJ32+8FGx52MZ8PsADl88zNhOY3ns7sdkOKS4NatXw759MH++2ZV4JEdC/hRQ3W67\nWuZr9k4C57XWyUCyUmod0AT4W8hHRET8tR4aGkpoaOitVSx8QmJyIh9aPuSHvT8wpM0QFnVbROHA\nwmaXJbxNejq8+y785z8Q7DvPBrBYLFgsFqd8liMXXgOBAxgXXuOAbUB3rXWM3T51gf8BYUAwsBXo\nqrXen+2z5MKrn0vPSGfKzikMtwync53OfNThIyoWr2h2WcJbff89TJpkjKjx4d8AXXrhVWudrpTq\nByzHNoQyRinV13hbT9JaxyqlfgeigHRgUvaAF8JyzEL/3/pTtkhZfuv1G00rNTW7JOHNrl2Df/8b\n5s3z6YAvKJlqWLjcHwl/8N6K99gRt4NPH/qUp+s9Lf3uouBGjoS9e2HOHLMrcTmXjpN3Jgl5/5J0\nI4nR60czccdEBrQcwIBWAygaVNTssoQviI+HBg2MB3TXrGl2NS7n6nHyQtySDJ3BjD0zGLp6KO1r\ntCfq9SiqlqpqdlnCl0REwPPP+0XAF5SEvHCqzSc20/+3/iilWPCPBbSs1tLskoSviYkxZpqUG58c\nIiEvnOLk5ZOErwxn7bG1jO44mp6NexKgZJJT4QLh4cbDucuVM7sSryD/CkWBXEu9xoi1I2g6sSk1\ny9Qktl8svZv0loAXrmGxQHS0MU+NcIicyYt80VozO3o24SvDaVmtJdv7bKdGmRpmlyV8WUYGDBwI\no0b51I1PriYhL25Z5KlI3v79bZLTkpn59Ez+747/M7sk4Q9mz4aAAOja1exKvIoMoRQOO3X5FENW\nDWHl0ZWM7DCSF5u+KN0ywj2Sk6FuXZg+Hdq2NbsatyvIEEr5FyryZO13bzyxMbeXup0D/Q7w8j0v\nS8AL9/nqK2jSxC8DvqCku0bkKkNnMGvvLIauGkqr21uxo88O6XcX7nfxIvz3v7B+vdmVeCUJeZGj\njX9uZMDyAWit+fGZH2ldvbXZJQl/NXIkPPus0V0jbpn0yYssjiYcZfDKwWw5uYWPO3ws492FuY4e\nhRYtjPniQ0LMrsY00icvCiwxOZFBKwbRfHJzGoc0lvHuwjO8/z707+/XAV9Q0l3j51LTU5m4fSIj\n14/kydpPEv1GNJVLVja7LCEgMhLWrYMpU8yuxKtJyPsprTU/H/iZ8JXh1ChTg5W9V9IopJHZZQlh\n0Breew8+/BCKFze7Gq8mIe+Htp7cysAVA7mUfIlxYePoVKuT2SUJkdUvv8D58/Dii2ZX4vUk5P3I\nkYtHGLp6KBv/3MiI9iN4ockLBAYEml2WEFmlpRmTkH3yCRSSiCoouarmB85dPUf/Zf25f8r9NK7Y\n+K+bmSTghUeaOhUqVYJHHzW7Ep8g/036sKs3rvL5ls/5YssXdG/Ynf1v7ZeHZgvPlpRkPBBkyRJ5\nbquTSMj7oNT0VKbsnMJH6z6iXY12bHl1C7XK1TK7LCHy9tln0L49NGtmdiU+Q0Leh2ToDOZEz2HY\nmmHULFuTJd2X0KyK/GMRXuLMGRg3DnbsMLsSnyJ3vPoArTW/HvqV91e/T3ChYEZ3HE2HOzuYXZYQ\nt+b116FECfj0U7Mr8TjyIG8/Zjlm4f3V75OYnMhH7T+iS90uKOnLFN7G+tzWAwfMrsTnSMh7qS0n\ntzBszTCOJhwlol0EPRr1kNEywnuFhxtNntvqdBLyXmbH6R18YPmAvfF7+Xfbf/NS05cICgwyuywh\n8m/tWti7F+bONbsSnyQh7yV2xu3kw7Ufsv30doa0GcJP//iJ4ELynEvh5azPbf34YyhSxOxqfJKE\nvIfbcXoHI9aNYPvp7YS3Dmf2M7MpGlTU7LKEcI65c415arp1M7sSnyWjazzU5hOb+WjdR0TFRzGo\n9SBeu/c1CXfhW1JSoF494w7X0FCzq/FoMrrGR2itWf3HakZtGMWRi0cY3GYwC7sulG4Z4ZvGj4cG\nDSTgXUzO5D1Ahs5g8YHFjN4wmkvJlwhvHU6vxr3kgqrwXQkJUKcOWCxQv77Z1Xi8gpzJS8ibKCUt\nhZl7Z/LJpk8oHlScIW2G8FTdp2QopPB9Awca89RMnGh2JV5BQt7LJFxPYNKOSYzbNo5GFRsxqPUg\n2tdoLzcxCf9w9Cg0b248t7VSJbOr8QrSJ+8lDl88zJdbvuSHvT/wRO0nWNpjKU0rNTW7LCHca+hQ\nePttCXg3kZB3Ma01q/5Yxbit49h0YhN9mvUh+o1oqpaqanZpQrjf1q2wfj18+63ZlfgN6a5xkcsp\nl/kh6gfGR44nQAXwrxb/omfjnhQLKmZ2aUKYQ2to2xZeegleftnsaryKdNd4kKj4KL6O/Jo5++bQ\n4c4OfPXIV4TWCJX+diEWLYJLl+CFF8yuxK84FPJKqTDgC4zHBX6rtf5vLvs1BzYBXbXWPzmtSg+X\ndCOJ2dGzmbxzMqcun+K1e19j7xt7pUtGCKsbN2DQIGNsfKCMHnOnPENeKRUAfAV0BE4DkUqpn7XW\nsTns9x/gd1cU6mm01mw8sZHvd3/PgpgFtL2jLcPaDiOsVhiFAuQXJCGy+PprqFULHn7Y7Er8jiNp\n1AI4pLU+DqCUmg10BmKz7fdPYD7Q3KkVepgjF48wc+9Mpu+ZTuHAwrzU9CX2vbmPKiWrmF2aEJ4p\nIcGYgGzNGrMr8UuOhHxV4ITd9kmM4P+LUqoK8JTWur1SKst7vuBM0hnm75/PzL0zOXLxCF0bdGXm\n0zNpUbWF9LULkZeRI6FLF2MKA+F2zupX+AIIt9v2+uSLuxLHwtiFzN03lz3xe3js7scY1nYYD9V8\nSKYbEMJRR47AtGnGjU/CFI6E/Cmgut12tczX7N0HzFbGae1twCNKqVSt9eLsHxYREfHXemhoKKEe\nMjmR1pqDFw7y84GfWRi7kNjzsTx696O80/IdOtXqRJFCMte1ELds8GB45x0ICTG7Eq9isViwWCxO\n+aw8x8krpQKBAxgXXuOAbUB3rXVMLvt/ByzJaXSNp42Tv556nfV/rmfpwaUsPbSU62nXebL2k3Sp\n14XQGqEUDixsdolCeK/166FXL4iNhaIyTXZBuHScvNY6XSnVD1iObQhljFKqr/G2npT9S/JTiDuk\nZaSx58weVv2xihVHV7Dl5BYahzTm0VqPMv8f82kS0kT62IVwhowM4wx+9GgJeJP59B2v11KvEXkq\nks0nN7Pu+Do2ndhE1VJVaV+jPQ/VfIjQGqGULlLabfUI4TemT4cJE2DzZpATpwKTWSiB5LRkos9G\nszNuJzvjdhJ5OpLY87E0qtiIltVa8n/V/4+2d7SlQvEKLvn+QohMV68ac8XPmwetWpldjU/wq5C/\nlHyJQxcPcfDCQQ5eOMi+c/uIPhvNHwl/ULt8bZpVaca9le41lpXvlQumQrhbRAQcPAizZpldic/w\n6pDXWnMj/QaXUi6RcD2BhOQELly7QPzVeM4kneFM0hlOXD7B8cTjHL90nJS0FO4ufzd3l7ub2uVr\n06BCAxpWbEjt8rXlMXlCmO3ECWjaFHbtgurV895fOMSrQv72sbeTrtNJy0jjeup1rqZeJVAFUiq4\nFGWLlqVc0XKUK1qOSiUqUal4JUJKhHB7qdu5o8wdVC9dnQrFKsjFUSE8Vc+ecNddMGKE2ZX4FK8K\n+WMJxyhFwg38AAAM10lEQVQUUIjAgECKBRWjaKGicnOREL5g0ybo2tUYMlm8uNnV+BSvCnlPGicv\nhHCSjAy4/37jiU89e5pdjc8pSMgHOLsYIYQfmj4dChWCHj3MrkRkI2fyQoiCuXLFGDK5aBG08Ln5\nCT2CdNcIIcwTHg5nzhgTkQmXkJAXQpjjwAFo3Rqio6FSJbOr8VnSJy+EcD+tjQutQ4ZIwHswCXkh\nRP4sWQLHjsE//2l2JeIm5GGkQohbl5xszDI5cSIUlim5PZmcyQshbt2nn0KTJvDQQ2ZXIvIgF16F\nELfm2DFo1gy2b4c77zS7Gr8gF16FEO7Tvz8MGCAB7yWkT14I4bglS4y5aebONbsS4SAJeSGEY65d\ng3/9CyZPhmCZ1ttbSHeNEMIxo0YZk5A9+KDZlYhbIBdehRB5s97ZumcPVK1qdjV+Ry68CiFcR2vo\n2xeGDZOA90IS8kKIm/v+e+Ph3P36mV2JyAfprhFC5O7cOWjYEH77De65x+xq/JbMQimEcI3evSEk\nxLjDVZimICEvQyiFEDlbuRLWr4d9+8yuRBSA9MkLIf4uKQn69IGvv5aHcns56a4RQvxd//6QkGA8\nu1WYTrprhBDOs3EjzJtnPO1JeD3prhFC2Fy/Di+/DP/7H5QrZ3Y1wgmku0YIYTN4MBw5YpzJC48h\n3TVCiILbsgW++w6iosyuRDiRdNcIIYwZJp9/Hr76yhgXL3yGdNcIIYyHcV+8CDNnml2JyIF01wgh\n8m/lSli0SLppfJR01wjhzxITjdE0U6ZA2bJmVyNcQLprhPBXWkOPHka4T5hgdjXiJlw+n7xSKkwp\nFauUOqiUCs/h/R5KqT2ZbYNSqlF+ihFCuNG0abB3L3z2mdmVCBfK80xeKRUAHAQ6AqeBSKCb1jrW\nbp+WQIzW+pJSKgyI0Fq3zOGz5ExeCE9w8KDxpKfVq6GRnJN5OlefybcADmmtj2utU4HZQGf7HbTW\nW7TWlzI3twDy+BghPNWNG0Y3zYcfSsD7AUdCvipwwm77JDcP8VeBZQUpSgjhQkOHQrVq8MYbZlci\n3MCpQyiVUu2Bl4A2ue0TERHx13poaCihoaHOLEEIcTOLFhlTFuzcCSpfv/0LN7BYLFgsFqd8liN9\n8i0x+tjDMrcHA1pr/d9s+zUGFgBhWusjuXyW9MkLYZYjR6BVK/jlF2jRwuxqxC1wdZ98JFBLKXWH\nUqow0A1YnK2A6hgB3zu3gBdCmOj6dXjmGRg+XALezzg0Tj5zxMyXGP8pfKu1/o9Sqi/GGf0kpdRk\n4GngOKCAVK31336S5ExeCJO8+ipcvQqzZkk3jReSB3kLIXI3cSKMGwfbtkGJEmZXI/JBQl4IkTOL\nBbp2NZ72VKuW2dWIfHL5Ha9CCC909Ch062Z00UjA+y0JeSF80ZUr0LkzvP8+dOxodjXCRNJdI4Sv\nSUszAr5qVfjmG7nQ6gOku0YIYdAa3nwT0tNh/HgJeCEPDRHCp3z8MWzfDmvXQlCQ2dUIDyAhL4Sv\nmDYNvv0WNm2CkiXNrkZ4CAl5IXzBokUQHg5r1kDlymZXIzyIhLwQ3u7336FPH1i2DOrVM7sa4WEk\n5IXwZmvXQq9e8PPP0KyZ2dUIDySja4TwVps2wXPPwezZ8MADZlcjPJSEvBDeaPVqYyz8jBlys5O4\nKQl5IbzN0qXGdAXz50OnTmZXIzychLwQ3mTePHj5ZViyBNq1M7sa4QXkwqsQ3uKLL+DTT43RNE2b\nml2N8BIS8kJ4uvR0ePddWLHCmDL4jjvMrkh4EQl5ITxZUhI8/zwkJBgBX6aM2RUJLyN98kJ4qkOH\njAdvlykDv/0mAS/yRUJeCE/0yy/QujW89ZYxH01wsNkVCS8l3TVCeJLUVBg+HKZPN+5ibdXK7IqE\nl5OQF8JTHD4MPXpAhQqwYweEhJhdkfAB0l0jhNm0NrpkWrWC3r2NrhoJeOEkciYvhJmOHIG+feHi\nRWOqgkaNzK5I+Bg5kxfCDKmp8MkncP/9EBYG27ZJwAuXkDN5Idxt2TIYMACqV4etW+Guu8yuSPgw\nCXkh3GXfPhg4EI4ehc8+g8cekwdtC5eT7hohXC021hg106EDPPww7N0Ljz8uAS/cQkJeCFeJijJG\ny7Rta/S3Hz4M77wDhQubXZnwIxLyQjiT1sYUBA89BI88AvXrG+E+ZAiULGl2dcIPSZ+8EM5w9ix8\n/z1MngzFihkXVrt3l7N2YToJeSHy6/p14ylNs2bBmjXQpYvxOL7775f+duExlNbafd9MKe3O7yeE\n0129CsuXw8KFxtOZmjUzztiffRZKlza7OuGjlFJorfN15iAhL8TNaG1M+btyJfz6K6xbZ5ypd+5s\nBHulSmZXKPyAhLwQzqI1HDwIGzbA+vXGVAMZGcaF1E6djIupcsYu3ExCXoj80BpOnYJdu2D7dqNt\n2wbFixtzubduDR07Qu3a0scuTOXykFdKhQFfYAy5/FZr/d8c9hkHPAJcBV7UWu/OYR8JeeF+GRlw\n8iQcOGC0mBjjhqToaAgKMh6K3bw53Hef0apVM7tiIbJwacgrpQKAg0BH4DQQCXTTWsfa7fMI0E9r\n/ZhS6n7gS611yxw+S0I+k8ViITQ01OwyPEKBj8W1axAXB6dPw4kT8OefxvKPP4x27BiULQt16kDd\nusayUSOjVazorD+GU8jPhY0cC5uChLwjQyhbAIe01sczv9lsoDMQa7dPZ2A6gNZ6q1KqtFIqRGsd\nn5+i/IH8ANv8dSy0NkavXL4MiYm2dvEiXLhga2fPGi0+Hs6cgeRk4wJolSrGpF+3324EeVgY1KwJ\nNWoYXTBeQH4ubORYOIcjIV8VOGG3fRIj+G+2z6nM1yTkPZ3WkJ5utLQ029LaUlNtS/t244ZtmZJi\nW1pbcrLRrl/P2q5dM4Lc2q5cMfrFv/zSWA8ONu4MLVPGOPsuW9ZYv+02KF8e6tWDdu2MM/CKFY2H\na5QrJ33mQuTC/TdDhYXd2v6Odu/cbD9H38u+X27vWdezL3N7L6fXTpwwhuTZ75Nby8jI/bWclhkZ\nRlhb1+1fs75uXdcaAgOztqAgKFQo63qhQsbdm0FBthYcbHstONjWihSxteBgI4iLFjVasWLGWbW1\nlSwJU6bAsGFQqpTxfYQQTuNIn3xLIEJrHZa5PRjQ9hdflVITgTVa6zmZ27FAu+zdNUop6ZAXQoh8\ncGWffCRQSyl1BxAHdAO6Z9tnMfAWMCfzP4XEnPrj81ukEEKI/Mkz5LXW6UqpfsBybEMoY5RSfY23\n9SSt9a9KqUeVUocxhlC+5NqyhRBCOMKtN0MJIYRwL5fMJ6+UClNKxSqlDiqlwnPZZ5xS6pBSardS\nqqkr6vAEeR0LpVQPpdSezLZBKeWzT3N25Ocic7/mSqlUpdTT7qzPnRz8NxKqlNqllIpWSq1xd43u\n4sC/kVJKqcWZWbFXKfWiCWW6nFLqW6VUvFIq6ib73Hpuaq2d2jD+4zgM3AEEAbuButn2eQRYmrl+\nP7DF2XV4QnPwWLQESmeuh/nzsbDbbxXwC/C02XWb+HNRGtgHVM3cvs3suk08FkOA0dbjAFwACpld\nuwuORRugKRCVy/v5yk1XnMn/dfOU1joVsN48ZS/LzVNAaaVUiAtqMVuex0JrvUVrfSlzcwvG/QW+\nyJGfC4B/AvOBs+4szs0cORY9gAVa61MAWuvzbq7RXRw5FhqwPlarJHBBa53mxhrdQmu9AUi4yS75\nyk1XhHxON09lD67cbp7yNY4cC3uvAstcWpF58jwWSqkqwFNa668BXx6J5cjPRW2gnFJqjVIqUinV\n223VuZcjx+IroL5S6jSwB+jvpto8Tb5yU+488RBKqfYYo5LamF2Lib4A7PtkfTno81IIuBfoABQH\nNiulNmutD5tblik6Abu01h2UUncBK5RSjbXWSWYX5g1cEfKngOp229UyX8u+z+157OMLHDkWKKUa\nA5OAMK31zX5d82aOHIv7gNlKKYXR9/qIUipVa73YTTW6iyPH4iRwXmudDCQrpdYBTTD6r32JI8fi\nJWA0gNb6iFLqD6AusN0tFXqOfOWmK7pr/rp5SilVGOPmqez/SBcDz8Nfd9TmePOUD8jzWCilqgML\ngN5a6yMm1OgueR4LrXXNzHYnRr/8mz4Y8ODYv5GfgTZKqUClVDGMC20xbq7THRw5FseBBwEy+6Br\nA0fdWqX7KHL/DTZfuen0M3ktN0/9xZFjAQwDygETMs9gU7XW2SeA83oOHossX+L2It3EwX8jsUqp\n34EoIB2YpLXeb2LZLuHgz8VI4Hu7oYWDtNYXTSrZZZRSs4BQoLxS6k9gOFCYAuam3AwlhBA+zCU3\nQwkhhPAMEvJCCOHDJOSFEMKHScgLIYQPk5AXQggfJiEvhBA+TEJeCCF8mIS8EEL4sP8HybcMa5B6\nLwwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a21a438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x_sing_plot, y_sing_plot_1[0])\n",
    "plt.plot(x_sing_plot, y_sing_plot_2[0])\n",
    "plt.plot(x_sing_plot, y_sing_plot_3[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 5. Solving with analytic Jacobians"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example 4 from [2]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Equations:\n",
    "$$\n",
    "y_1' = 3 T (y_1 + y_2 - \\frac{1}{3} y_1^3 - 1.3) \\\\\n",
    "y_2' = -T (y_1 - 0.7 + 0.8 y_2) / 3\n",
    "$$\n",
    "Boundary conditions:\n",
    "$$\n",
    "y_1(0) - y_1(1) = 0 \\\\\n",
    "y_2(0) - y_2(1) = 0 \\\\\n",
    "T (y_1(0) - 0.7 + 0.8 y_2(0)) + 3 = 0\n",
    "$$\n",
    "$T$ is an unknown parameter (period in the original variables)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's write Jacobian matrices:\n",
    "$$\n",
    "\\frac{\\partial f}{\\partial y} = \\begin{pmatrix}\n",
    "3 T (1 - y_1^2) & 3 T \\\\\n",
    "-T / 3 & -0.8 T / 3\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "$$\n",
    "\\frac{\\partial f}{\\partial T} = \\begin{pmatrix}\n",
    "3 (y_1 + y_2 - \\frac{1}{3} y_1^3 - 1.3) \\\\\n",
    "-(y_1 - 0.7 + 0.8 y_2) / 3\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "$$\n",
    "\\frac{\\partial bc}{\\partial y_a} = \\begin{pmatrix}\n",
    "1 & 0 \\\\\n",
    "0 & 1 \\\\\n",
    "T & 0.8 T\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "$$\n",
    "\\frac{\\partial bc}{\\partial y_b} = \\begin{pmatrix}\n",
    "-1 & 0 \\\\\n",
    "0 & -1 \\\\\n",
    "0 & 0\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "$$\n",
    "\\frac{\\partial bc}{\\partial T} = \\begin{pmatrix}\n",
    "0 \\\\ 0 \\\\ y_1(0) - 0.7 + 0.8 y_2(0)\n",
    "\\end{pmatrix}\n",
    "$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fun_imp(x, y, p):\n",
    "    T = p[0]\n",
    "    return np.vstack((\n",
    "        3 * T * (y[0] + y[1] - y[0]**3/3 - 1.3),\n",
    "        -T / 3 * (y[0] - 0.7 + 0.8 * y[1])\n",
    "    ))\n",
    "\n",
    "def bc_imp(ya, yb, p):\n",
    "    T = p[0]\n",
    "    return np.array([ya[0] - yb[0], \n",
    "                     ya[1] - yb[1],\n",
    "                     T * (ya[0] - 0.7 + 0.8 * ya[1]) + 3])\n",
    "\n",
    "\n",
    "def fun_jac_imp(x, y, p):\n",
    "    T = p[0]\n",
    "    df_dy = np.empty((2, 2, x.shape[0]))\n",
    "    df_dy[0, 0] = 3 * T * (1 - y[0]**2)\n",
    "    df_dy[0, 1] = 3 * T\n",
    "    df_dy[1, 0] = -T / 3\n",
    "    df_dy[1, 1] = -0.8 * T / 3\n",
    "    \n",
    "    df_dp = np.empty((2, 1, x.shape[0]))\n",
    "    df_dp[0] = 3 * (y[0] + y[1] - y[0]**3/3 - 1.3)\n",
    "    df_dp[1] = -(y[0] - 0.7 + 0.8 * y[1]) / 3\n",
    "    \n",
    "    return df_dy, df_dp\n",
    "\n",
    "def bc_jac_imp(ya, yb, p):\n",
    "    T = p[0]\n",
    "    dbc_dya = np.array([\n",
    "        [1, 0],\n",
    "        [0, 1],\n",
    "        [T, 0.8 * T]\n",
    "    ])\n",
    "    dbc_dyb = np.array([\n",
    "        [-1, 0],\n",
    "        [0, -1],\n",
    "        [0, 0]\n",
    "    ])\n",
    "    dbc_dp = np.array([\n",
    "        [0],\n",
    "        [0],\n",
    "        [ya[0] - 0.7 + 0.8 * ya[1]],\n",
    "    ])\n",
    "    \n",
    "    return dbc_dya, dbc_dyb, dbc_dp\n",
    "\n",
    "x_imp = np.linspace(0, 1, 5)\n",
    "\n",
    "y_imp = np.empty((2, 5))\n",
    "y_imp[0] = np.sin(2 * np.pi * x_imp)\n",
    "y_imp[1] = np.cos(2 * np.pi * x_imp)\n",
    "\n",
    "p_imp = np.array([2 * np.pi])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Iteration    Max residual    Total nodes    Nodes added  \n",
      "       1          4.16e+00           5              6       \n",
      "       2          2.68e+00          11             12       \n",
      "       3          2.89e-01          23              8       \n",
      "       4          7.42e-02          31              5       \n",
      "       5          9.53e-03          36              0       \n",
      "Solved in 5 iterations, number of nodes 36, maximum relative residual 9.53e-03.\n"
     ]
    }
   ],
   "source": [
    "res_imp = solve_bvp(fun_imp, bc_imp, x_imp, y_imp, p=p_imp, fun_jac=fun_jac_imp,\n",
    "                    bc_jac=bc_jac_imp, tol=1e-2, verbose=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T=10.71 --- matches with the value from the paper\n"
     ]
    }
   ],
   "source": [
    "print(\"T={:.2f} --- matches with the value from the paper\".format(res_imp.p[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_plot_imp = np.linspace(0, 1, 100)\n",
    "y_plot_imp = res_imp.sol(x_plot_imp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10a851ba8>]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm0VNW17/HvpLcltqgQAUFETZRrgxiIlIiK3kQlGgWM\nGtRAkqHEZzRqjAbHTW4CZhhjk4coGGwQjU1EsQGb0pgoIp0cehQUUOEhAQVFDzDfH+ugiOdwmtpV\nq2rX7zNGDapOLfae7HHYc6/e3B0RESlPjWIHICIi8SgJiIiUMSUBEZEypiQgIlLGlARERMqYkoCI\nSBnLOQmYWRsze8HMZpvZLDMbUkO5W8xsoZnNMLMuuZ5XRERy1ySBY2wELnf3GWa2MzDVzCa6+7wt\nBczsFKCDux9oZscAI4BuCZxbRERykHNNwN0/cPcZVe/XAXOB1tsUOx24p6rMZKClmbXK9dwiIpKb\nRPsEzKwd0AWYvM1XrYGlW31eztcThYiIFFhiSaCqKehh4BdVNQIRESlySfQJYGZNCAngXnd/vJoi\ny4FvbvW5TdXPqjuWFjMSEaknd7eG/L2kagKjgTnu/pcavh8PnA9gZt2ANe6+oqaDubte7vz2t7+N\nHkMxvHQddC10Lbb/ykXONQEz6w6cC8wys+mAA78G2ob7uY9096fM7FQzWwSsBwbmel4REcldzknA\n3f8FNK5DuUtyPZeIiCRLM4aLWCaTiR1CUdB1+JKuxZd0LZJhubYnJc3MvNhiEhEpZmaGR+4YFhGR\nEqQkICJSxpQERETKmJKAiEgZUxIQESljiSwbIeVnzRp44gnYvBnat4cDDoD99oNGeqwQKSlKAlJn\n7vDYY3DvvfDCC9CrF+y0EyxeHF6bNsE558B558FRR4E1aMCaiBSS5glInd12W3hdfTX07QstW371\n+7fegvvuC0miaVP4P/8HLrgAmjePE69IuchlnoCSgNTJ5Mnw/e/Dq69Chw7bL+sOL78Mw4bBzJlw\n+eUweDDsvHNhYhUpN5osJnn14YehmWfkyNoTAIRmoJ494amn4MknQwI58MBQi/j88/zHKyJ1p5qA\nbNfmzaEGcPDB8Kc/Nfw4M2bANdfA/Pnwu99Bv37qRBZJipqDJG9GjAht/NlsaOfPVTYLv/pVeP/n\nP0P37rkfU6TcqTlI8mbUKLjhhmQSAEAmA6+9BkOGQP/+cPbZsGRJMscWkfpTEpAaLVgAy5bB8ccn\ne9xGjeBHP4J58+Bb34IjjwyJ5tNPkz2PiNROSUBq9MAD4Um9ca1bBjXMjjvC9dfDtGkwe3bod3j0\n0TC6SEQKQ30CUi136NwZ7rkHjjmmMOd84QW49FLYf3+45ZYwokhEaqc+AUnc9OmwcSN07Vq4c/bq\nFUYRnXACHHssXHedmohE8i2RJGBmo8xshZm9WcP3Pc1sjZlNq3r9JonzSv6MHQsDBhR+6YemTeGK\nK0IymDcPvv1tmDixsDGIlJNEmoPMrAewDrjH3Q+r5vuewC/d/bQ6HEvNQZFt3hyaZCZNCu30MT31\nFFxySaiR3Hwz7LNP3HhEilH05iB3fwX4Ty3FtJxYifjnP2HPPeMnAIBTT4WKirBS6WGHwZ13hiQl\nIskoZJ/AsWY2w8wmmNkhBTyv1NOWpqBiseOO8Ic/wHPPwV13hSGr8+fHjkokHQqVBKYC+7t7F+A2\n4B8FOq/UkzuMHw9nnRU7kq877DD4979DbD16hMRQWRk7KpHSVpD9BNx93Vbvnzazv5rZ7u6+urry\nQ4cO/eJ9JpMhk8nkPUYJli0LiaB9+9iRVK9x4zCM9LTTYNAg+PvfYfRo6NIldmQihZPNZslms4kc\nK7F5AmbWDnjC3b9dzXet3H1F1fuuwEPu3q6G46hjOKJHH4W77w67hhU7dxgzJqxF9LOfwbXXQrNm\nsaMSKbzoHcNmNhb4N9DJzN41s4FmNtjMBlUVOcvMKsxsOnAzcE4S55Xkvf56YecG5MIMfvzjMJx0\n2rQQ98yZsaMSKS2aMSxf0atXeLLu0yd2JPXjHmY3X3llWJzu6quhiTZPlTKhpaQlEZs3w267wdtv\nwx57xI6mYZYtg4EDYd26kBS09ISUg+jNQZIO8+fDXnuVbgIAaNMGnn02DHH9znfgjju0IJ3I9igJ\nyBdefx2OPjp2FLlr1CiMIHr55ZAE+vaFVatiRyVSnJQE5Aul1ClcFwcfDK++Ch07hiGkL7wQOyKR\n4qMkIF+YMiVdSQCgefOwN/KoUWEjm9/8JqyOKiKBOoYFgM8+C53Cq1aFZRrSaMUKOPfckATGjoX9\n9osdkUgy1DEsOZs5Ezp1Sm8CAGjVKnQa9+oVtrScNCl2RCLxKQkIEJqC0tApXJvGjcOWlvffDxdc\nAH/8o0YPSXlTEhAgfZ3CtenVK/yb//EPOPNM+Oij2BGJxKEkIED5JQEIcwpeegn23jv82xcujB2R\nSOEpCQhr18LSpXDoobEjKbzmzWHECLj88rA8tYaRSrlREhCmTw9r9ZfzWjuDBsG4cWGm8YgRsaMR\nKRwlAaGiIiSBcnf88fDKK/CXv4SagbaxlHKgJCBUVJRnU1B1OnYMu5e98Qb07x/mT4ikmZKAMHs2\nfOtbsaMoHrvtBhMnwqZNYUnttWtjRySSP0oCZc5dNYHqtGgBDz4YkuNxx8HKlbEjEskPJYEy9/77\n0LRpGCYpX9W4MdxyC5xxBvTsCcuXx45IJHllPB5EQLWA2pjBDTeE5TR69oTnn4e2bWNHJZIcJYEy\np/6AurnqqpAIjjsOXnwRDjggdkQiyUhqo/lRZrbCzN7cTplbzGyhmc0wsy5JnFdyN3u2agJ1deml\nYe/i3r3DNpYiaZBUn8DdwMk1fWlmpwAd3P1AYDCg6ThFoqJCNYH6+NnP4Oc/D4lgxYrY0YjkLpEk\n4O6vAP/ZTpHTgXuqyk4GWppZqyTOLQ3nDnPmqCZQX1dcEeYQnHQSrF4dOxqR3BRqdFBrYOlWn5dX\n/Uwievdd2HnnMC5e6uf66+HEE+H734cNG2JHI9JwRdkxPHTo0C/eZzIZMplMtFjSTJ3CDWcGw4eH\nGsHAgWF/gkYacC0Fks1myWaziRwrse0lzawt8IS7f20VGjMbAbzo7g9WfZ4H9HT3r7WqanvJwhk+\nPMwT+POfY0dSuj79NOxN0Ls3/M//xI5GylWxbC9pVa/qjAfOBzCzbsCa6hKAFJZqArnbYQd4/PFQ\nExgzJnY0IvWXSHOQmY0FMsAeZvYu8FugGeDuPtLdnzKzU81sEbAeGJjEeSU3FRVhtIvkZu+9YcKE\nMJnskEPKY5tOSY/EmoOSouagwti8GXbZJTQH7bpr7GjS4ZFHwsihadPU2S6FVSzNQVJCFi+GPfdU\nAkjSmWfC6aeHDez1HCOlQkmgTGmSWH4MHx5WHP3Tn2JHIlI3RTlEVPJPy0XkR7Nm8NBDYeP6Hj3g\n2GNjRySyfaoJlCmtHpo/++8Pt98e5g98+mnsaES2T0mgTM2aBd/+duwo0uvMM8P13Wreo0hR0uig\nMvT559CyZVj3ZocdYkeTXitXwmGHwfjxoXlIJF80OkjqZcGCsDGKEkB+7b13mI194YXasF6Kl5JA\nGVJTUOH06wcdO8Lvfx87EpHqKQmUISWBwjELncS33x5WbRUpNkoCZWjWLM0RKKTWrcPyHNddFzsS\nka9Tx3AZat8eJk6EAw+MHUn5+Ogj6NQJnnkGumhzVUlYLh3DSgJl5uOPYd99Ye1aaNw4djTl5fbb\nw4qjEyfGjkTSRqODpM4qKuDgg5UAYhg0CJYsgWefjR2JyJeUBMqM+gPiadoU/vhHuPJK2LQpdjQi\ngZJAmamo0MigmPr2hebN4cknY0ciEigJlBkND43LDH75S7jpptiRiARKAmXEXUmgGJx1VugbeOON\n2JGIKAmUlQ8+CH+2ahU3jnLXpAkMGaLagBQHJYEysqU/wBo0kEySdPHFYc6AZhFLbIkkATPrY2bz\nzGyBmV1Vzfc9zWyNmU2rev0mifNK/agpqHi0bAk//jHcemvsSKTc5ZwEzKwRcBtwMnAo0N/MOldT\n9GV3P6Lq9btczyv1pyRQXIYMgdGjwwQ+kViSqAl0BRa6+zvuXgmMA06vppwaISJTEigu7dpB797w\nt7/FjkTKWRJJoDWwdKvPy6p+tq1jzWyGmU0ws0MSOK/Uw6ZNMHeutpQsNhdeCPfdFzsKKWeF2mh+\nKrC/u39iZqcA/wA61VR46FZ78mUyGTKZTL7jS705c6BNG9hll9iRyNZOOAHOPx8WLQr7DojURTab\nJZvNJnKsnBeQM7NuwFB371P1+WrA3X3Ydv7OYuBId19dzXdaQC4PRo+GF17QU2cxuvTSMGz3Nxou\nIQ0UewG5KUBHM2trZs2AfsD4bQJstdX7roTk87UEIPkzZYr2uS1WAwbA/feHyXwihZZzEnD3TcAl\nwERgNjDO3eea2WAzG1RV7CwzqzCz6cDNwDm5nlfq5/XX4eijY0ch1enWDTZsgJkzY0ci5Uj7CZSB\nDRtg993hww+1uXyxuvZaqKyE4cNjRyKlKHZzkBS5mTOhc2clgGI2YAA88ABs3hw7Eik3SgJlQE1B\nxe/QQ0Nt7ZVXYkci5UZJoAxMmaIkUAoGDICxY2NHIeVGSaAMaGRQaTj7bHjsMTUJSWEpCaTc2rWw\ndCkcojnaRa99e9h117Daq0ihKAmk3NSp0KVLWMNeil/v3vDcc7GjkHKiJJBy6hQuLUoCUmhKAimn\n/oDScvzxYYTQ55/HjkTKhZJAymlkUGnZfXc46CB47bXYkUi5UBJIsQ8+gHXroEOH2JFIfahJSApJ\nSSDFttQCtKdwaVESkEJSEkixbBZ69IgdhdRX9+5hF7i1a2NHIuVASSDFJk2CE0+MHYXUV4sWYWXR\nl16KHYmUAyWBlPrggzBJ7KijYkciDaEmISkUJYGUeu65MNxQk8RKk5KAFIqSQEpNmgQnnRQ7Cmmo\nLl1gxQpYvjx2JJJ2SgIp5K7+gFLXuHHoINZ8Ack3JYEUmj07dC5qfkBpO/JImDYtdhSSdkoCKaRa\nQDoccYSSgORfIknAzPqY2TwzW2BmV9VQ5hYzW2hmM8ysSxLnleopCaTDEUeEVWC15bbkU85JwMwa\nAbcBJwOHAv3NrPM2ZU4BOrj7gcBgYESu55XqffZZWICsV6/YkUiu9tsvzPZ+773YkUiaJVET6Aos\ndPd33L0SGAecvk2Z04F7ANx9MtDSzFolcG7ZxquvwsEHh4XIpLSZfVkbEMmXJJJAa2DpVp+XVf1s\ne2WWV1NGEqCmoHRRv4DkW1FOJRo6dOgX7zOZDJlMJlospWbCBLj11thRSFKOOALGjIkdhRSbbDZL\nNptN5FjmOfY6mVk3YKi796n6fDXg7j5sqzIjgBfd/cGqz/OAnu6+oprjea4xlas5c8IEsXffhUYa\n95UKixfDd78Ly5bFjkSKmZnh7g1aLziJW8UUoKOZtTWzZkA/YPw2ZcYD58MXSWNNdQlAcnP//dC/\nvxJAmrRrB598EmYPi+RDzrcLd98EXAJMBGYD49x9rpkNNrNBVWWeAhab2SLgDuDnuZ5Xvsodxo6F\nc8+NHYkkaUvn8PTpsSORtEqkT8DdnwEO2uZnd2zz+ZIkziXV+/e/Yccd4fDDY0ciSdvSOdynT+xI\nJI3UcJAS998fagHaRSx9NExU8innjuGkqWO4/iorw8SiKVNCG7Kky4IFcPLJoZNYpDqxO4Ylsmef\nhYMOUgJIq44d4cMPYfXq2JFIGikJpMCWpiBJp0aNwv4C6hyWfFASKHEffwxPPw0//GHsSCSftKy0\n5IuSQIm7664wQWzPPWNHIvl0+OEwc2bsKCSNlARK2Oefw003wVXVLt4tadKpEyxcGDsKSSMlgRL2\nwAOhQ/jII2NHIvnWqVMYJaSBc5I0JYEStXkzDB+uWkC52GOPMAdk1arYkUjaKAmUqAkToHlz6N07\ndiRSCGZw4IFqEpLkKQmUqGHDQi1AM4TLx5YmIZEkKQmUoH/9K2w5eOaZsSORQlISkHxQEigx7nDt\ntXD11dCkKLcEknxRc5Dkg5JAifn732HNGrjootiRSKGpJiD5oAXkSsj69WET+fvug+OOix2NFNrH\nH8M++4Q/tXGQbE0LyJWJYcOge3clgHK1yy6w666hP0gkKWpVLhGLF8Ptt8OMGbEjkZi2NAm1aRM7\nEkkL1QRKgDtcfnl4ffObsaORmNQvIElTTaAE3HsvzJ8flomQ8qYRQpK0nJKAme0GPAi0BZYAZ7v7\n2mrKLQHWApuBSnfvmst5y8miRfDLX8Lzz0OLFrGjkdg6dYJ//jN2FJImuTYHXQ085+4HAS8A19RQ\nbjOQcff/UgKou88/h/794frr4bDDYkcjxUDNQZK0nIaImtk8oKe7rzCzfYCsu3euptxi4Ch3/7AO\nx9QQ0SpXXQVz5sD48VoeQoING+Ab34B16zRZUL4Uc4jo3u6+AsDdPwD2rqGcA5PMbIqZ/STHc5aF\nCRPCfIDRo5UA5EstWsC++8KSJbEjkbSo9VnCzCYBrbb+EeGm/ptqitf0CN/d3d83s70IyWCuu79S\n0zmHDh36xftMJkMmk6ktzFSZOhUGDgw1gL32ih2NFJstTUIdO8aORGLJZrNks9lEjpVrc9BcQlv/\nluagF9394Fr+zm+Bj939phq+L+vmoCVLwoSw22+HM86IHY0Uo0suCQngsstiRyLFImZz0Hjgx1Xv\nLwAe37aAme1oZjtXvd8JOAmoyPG8qbR6NZxySlgcTglAaqKtJiVJuSaBYcCJZjYfOAH4I4CZ7Wtm\nT1aVaQW8YmbTgdeAJ9x9Yo7nTZ2PPoLvfQ9OPRUuvTR2NFLMNEJIkqQF5IrAmjVw8slw1FFw661a\nHEy27623oFcveOed2JFIscilOUhJILLVq+Gkk6BHD/jznzUSSGq3cSPstFOoPTZvHjsaKQZaRbRE\nffABnHACHH+8EoDUXZMm0Lo1LF0aOxJJAyWBSCoqoFu3sEXk8OFKAFI/7dpproAkQ3MOI5g0Cc49\nF26+GQYMiB2NlCIlAUmKagIF5B7G/593HjzyiBKANFzbtkoCkgzVBArk449h0CCYOxdeeUWzPSU3\n7dqFGqVIrlQTKICKCjj6aNh5Z3j1VSUAyZ2agyQpSgJ5tHkz3HQTZDLw61/DnXfCDjvEjkrSoF07\nzROQZKg5KE/efjssArd5M0yeDB06xI5I0qR1a1i5Muw50axZ7GiklKkmkLCNG8Oon2OOgdNPh2xW\nCUCS16QJ7Lef5gpI7lQTSNDkyfDTn8Luu4fO34MOih2RpNmWfgE9ZEguVBNIwMqVYeTPGWfAFVfA\nc88pAUj+qXNYkqAkkIPPPoMbb4RDDglrucyZEyaBafavFILmCkgSlAQawB0eeggOPRRefhn+9a+w\n9s9uu8WOTMqJagKSBPUJ1NOLL8KvfhUSwR13hAXgRGLQMFFJgpJAHb3+Olx3HSxaBL//PZx9ttb9\nl7hUE5Ak6DZWi5kzw1DPM8+EH/wgLPvQr58SgMTXpg2sWBHmCog0lG5lNZg6NYz26dMnrPe/cCEM\nHqyJOVI8mjSBffeFZctiRyKlTElgG6++Gvb6Pe20sIXf22/DZZdBixaxIxP5OjUJSa5ySgJmdpaZ\nVZjZJjM7Yjvl+pjZPDNbYGZX5XLOfHCHZ56Bnj3D8s6nnhr2cR0yRGv9SHHTMFHJVa4dw7OAvsAd\nNRUws0bAbcAJwHvAFDN73N3n5XjunFVWhqGef/pTeH/11aG9v4m6y6VEqCYgucrpdufu8wHMtjs9\nqiuw0N3fqSo7DjgdiJYE1q6Fu+6Cv/wFDjwQ/vd/Q9u/JnlJqWnXLqxPJdJQhegTaA1svczVsqqf\nFdxbb8EvfgHt24eO38ceg+efh1NOUQKQ0qSagOSq1pqAmU0CWm39I8CBa939iXwENXTo0C/eZzIZ\nMplMg4/lHiZ43XJLWNTt4ovhzTfD8DqRUqckUJ6y2SzZhKqA5u65H8TsReCX7j6tmu+6AUPdvU/V\n56sBd/dhNRzLk4hp/Xq47z649daQCIYMgR/9KKzxI5IWlZVhx7p166Bp09jRSCxmhrs3qD0jyeag\nmgKYAnQ0s7Zm1gzoB4xP8LxfsWgRXH457L8/PP10WNu/oiKM8VcCkLRp2hT22UdzBaThch0ieoaZ\nLQW6AU+a2dNVP9/XzJ4EcPdNwCXARGA2MM7d5+YW9ldt2gRPPBE6d489NkzomjoV/vEP6N1b7f2S\nbhomKrlIpDkoSfVpDlqxAkaNCgu57bsv/Pzn8MMfamy/lJfzzw/7WF94YexIJJZcmoNKbkS8exgS\nN2IETJwIZ50Fjz4KRx4ZOzKRODp0CDPbRRqiZJLAqlUwZgyMHBnaQX/60/C+ZcvYkYnE1aEDTJgQ\nOwopVUWdBNzDpi0jR8JTT4X1fO6+O7T7q51fJOjQIcyBEWmIokwCK1fCPffAnXeGp/5Bg8JQz913\njx2ZSPFREpBcFGXHcMuWTt++YWLXd76jp36R7XGHXXeFpUvhG9+IHY0U2oYNsMMOKesYfucdtfWL\n1JXZl7UBDZAoH2++GdZAGzs2t+MU5X4CSgAi9aMmofLw0Uehj7RrV/jv/4bddoM33sjtmEVZExCR\n+lESSC/3sO7ZqFFhAuwJJ8ANN8BJJ0HjxrkfX0lAJAUOOCD3J0IpLu+/HwbIjB4d9jS/6CIYPhz2\n3jvZ8xRlc5CI1I9qAulQWRme9k87DQ45JOxt/re/wZw5cMUVyScAUE1AJBWUBErb3Lnhif/ee6Fj\nx/DUP3ZsWCE235QERFJg//3DWlqffQbNm8eORurio4/gwQfDzf+dd+CCC8Lk2E6dChuHkoBICjRp\nAt/8JixeDJ07x45GarJlFYTRo+Hxx6FXL7j22rACcqy9zZUERFJiS5OQkkDxWbYsrH12993QokVY\n8fXGG/PTxl9fSgIiKaF+geKyYQOMHx+e+l9/Hc45Bx54AI46qrhWQVASEEkJJYHiMH16uPE/8AAc\nfnh46n/0Udhxx9iRVU9JQCQlOnSAF1+MHUV5+vDDMJpn9GhYvRoGDoQpU6B9+9iR1U5JQCQlVBMo\nrE2bYNKkcOOfODEs43DjjaGzt1EJzcDKaRVRMzsLGAocDBzt7tNqKLcEWAtsBirdvet2jlnn7SVF\n5Evr18Oee4Y/S+kmVGoWLQodvGPGwH77heaefv3iruAac3vJWUBf4I5aym0GMu7+nxzPJyI12Gmn\ncCNavjwMF5XkrF8PDz8cnvrnzoXzzoNnnoFvfSt2ZLnLKQm4+3wAs1r7ug0tUSGSd1uahJQEcucO\nkyeHhdsefhi6d4fLLgvNPs2axY4uOYXqE3BgkpltAka6+50FOq9IWdmSBDKZ2JGUrhUrwvINo0fD\nxo2hk3f27ND0k0a1JgEzmwS02vpHhJv6te7+RB3P093d3zezvQjJYK67v1L/cEVke9Q53DAbN8LT\nT4en/pdegjPOCOv2d+9eXGP686HWJODuJ+Z6End/v+rP/2dmjwFdgRqTwNChQ794n8lkyOixRqRO\nOnQIE5SkbubPD0/899wThnNedFGoBeyyS+zIti+bzZLNZhM5ViJ7DJvZi8AV7j61mu92BBq5+zoz\n2wmYCNzg7hNrOJZGB4k00LRpodNy9uzYkRSv9evh738PT/0LF4brdeGFcPDBsSNruFxGB+U6RPQM\n4FZgT2ANMMPdTzGzfYE73f17ZtYeeIzQhNQEuN/d/7idYyoJiDRQZWXYcvC998Lm8xK4h6Ubtu7k\nveii0MnbtGns6HIXLQnkg5KASG6++10YOjRsQ1juVq2C++4LN/8NG8IT/wUXpK+TN5ckoGGbIilz\nzDHw2muxo4hn8+Ywg/ecc8IGLVOnwm23wYIFcM016UsAudKyESIp061bmM1abpYuDTN5R4+GPfYI\nzT133BF3Jm8pUHOQSMosWwb/9V+wcmX6hzdWVsITT8Bdd4WJXf36hZv/EUfEjqywYi4bISJFpk2b\nsMXk22+HIaNptHBhuPGPGQMHHRRu/A8/XLzLNRcz9QmIpFC3bunrF9iwAe6/P8yG7tEjtP2/9FJ4\nnX++EkBDKQmIpFC3bqF5JA0qKuAXvwg1nHvugUsuCe3/N94YagGSGzUHiaRQt27w0EOxo2i4Tz4J\n8Y8cCe+8E4Z2vvEGtGsXO7L0UcewSAp98gnstVcYJ7/DDrGjqbuZM8ON/4EHwoSun/wETj0Vmuhx\ndbvUMSwiX7HjjtC5c9jv9jvfiR3N9q1fD+PGhZv/e+/BxReHZKDlsAtDSUAkpbZ0DhdrEpg5M4zj\nHzcuzHK+/nro0wcaN44dWXlREhBJqW7dwhj6YrJ+PTz4YLj5v/9+eOp/883Q6StxqE9AJKUWLoTe\nvUPHamwVFeHGP3ZsqJkMHgynnKKn/qSoT0BEvqZjx/Dk/fbbcMABhT//hg1hAteIEbB4cXjqnz4d\n9t+/8LFIzVQTEEmxX/8aVq8ON+JCWbAgdPKOGQNHHgk//Sl873sa4ZNPWkpaRKq1alWYUDVtGrRt\nm7/zVFbC44+HZDNrVtiX9yc/Se+yFcVGSUBEapTP2sCSJXDnnWHlzk6d4Gc/g759w9pFUjhKAiJS\no6RrAxs3woQJoaP39dfhRz8KHb2lvD1jqVMSEJHtSqI2sGRJWK9/1KjQuTt4MJx9dmnNSE4rJQER\n2a6G1gY+/RTGjw83/mnTYMCAMMrnsMPyF6vUX7TtJc1suJnNNbMZZvaImVW7tbWZ9TGzeWa2wMyu\nyuWcIlJ/e+4JQ4aEGbmTJm2/7KZNkM2GNfpbtw5t/gMHhs1qbrlFCSBtcl1KeiJwqLt3ARYC12xb\nwMwaAbcBJwOHAv3NrHOO5y0L2Ww2dghFQdfhS7lci+uvh2HDQuftaafBjBnw1lswezZMmQJ//Suc\neWZYeO6yy0Ib/6xZ8Nxz0L8/tGiR3L8jCfq9SEZOScDdn3P3zVUfXwOqm/zdFVjo7u+4eyUwDjg9\nl/OWC/2SB7oOX8rlWpiFm//s2WGtnr594cQT4Yc/hEGDQidv375hdu+MGXDFFaEmUKz0e5GMJKdv\nXEi4wW+aCbRfAAAD20lEQVSrNbB0q8/LCIlBRCJo3hyuvDK8RGpNAmY2CWi19Y8AB6519yeqylwL\nVLr72LxEKSIieZHz6CAz+zHwE6CXu39WzffdgKHu3qfq89WAu/uwGo6noUEiIvUUZQE5M+sDXAkc\nV10CqDIF6GhmbYH3gX5A/5qO2dB/iIiI1F+uo4NuBXYGJpnZNDP7K4CZ7WtmTwK4+ybgEsJIotnA\nOHefm+N5RUQkAUU3WUxERAon15pAg9Rl8piZ3WJmC6smonUpdIyFUtu1MLMBZjaz6vWKmX07RpyF\nUNdJhWZ2tJlVmtkPChlfIdXx/0jGzKabWYWZvVjoGAulDv9HdjWz8VX3illV/ZSpZGajzGyFmb25\nnTL1u3e6e0FfhMSzCGgLNAVmAJ23KXMKMKHq/THAa4WOs4iuRTegZdX7PuV8LbYq9zzwJPCD2HFH\n/L1oSWhebV31ec/YcUe8FtcAf9hyHYAPgSaxY8/T9egBdAHerOH7et87Y9QE6jJ57HTgHgB3nwy0\nNLNWpE+t18LdX3P3tVUfXyPMu0ijuk4qvBR4GFhZyOAKrC7XYgDwiLsvB3D3VQWOsVDqci0c2KXq\n/S7Ah+6+sYAxFoy7vwL8ZztF6n3vjJEEqps8tu2Nbdsyy6spkwZ1uRZbuxh4Oq8RxVPrtTCz/YAz\n3P3/EuarpFVdfi86Abub2YtmNsXMzitYdIVVl2txG3CImb0HzAR+UaDYilG9753a8K1EmNnxwEBC\ndbBc3Qxs3Sac5kRQmybAEUAvYCfgVTN71d0XxQ0ripOB6e7ey8w6EEYrHubu62IHVgpiJIHlwNZb\nTbep+tm2Zb5ZS5k0qMu1wMwOA0YCfdx9e1XBUlaXa3EUMM7MjND2e4qZVbr7+ALFWCh1uRbLgFXu\nvgHYYGYvA4cT2s/TpC7XYiDwBwB3f8vMFgOdgTcKEmFxqfe9M0Zz0BeTx8ysGWHy2Lb/iccD58MX\nM47XuPuKwoZZELVeCzPbH3gEOM/d34oQY6HUei3c/YCqV3tCv8DPU5gAoG7/Rx4HephZYzPbkdAJ\nmMb5N3W5Fu8AvQGq2r87AW8XNMrCMmquBdf73lnwmoC7bzKzLZPHGgGj3H2umQ0OX/tId3/KzE41\ns0XAekKmT526XAvgOmB34K9VT8CV7p66BfjqeC2+8lcKHmSB1PH/yDwzexZ4E9gEjHT3ORHDzos6\n/l78DvjbVsMmf+XuqyOFnFdmNhbIAHuY2bvAb4Fm5HDv1GQxEZEyFmWymIiIFAclARGRMqYkICJS\nxpQERETKmJKAiEgZUxIQESljSgIiImVMSUBEpIz9fz4aIHjN2SwYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x104b799e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x_plot_imp, y_plot_imp[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the starting guess (sinusoid) was poor and the problem turned out to be hard. If you set 10 nodes in the initial mesh, the solver won't converge to the correct solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}