{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nonlinear Equations\n",
    "\n",
    "## Overview, Objectives, and Key Terms\n",
    "\n",
    "In [Lesson 19](L19_NumericalDifferentiation.ipynb), finite differences gave us practical derivative estimates, including central differences and complex-step derivatives for callable functions. Here, derivatives become part of an iteration: we solve nonlinear equations by repeatedly asking how the residual changes.\n",
    "\n",
    "A scalar nonlinear equation has the form\n",
    "\n",
    "$$\n",
    "  f(x)=0.\n",
    "$$\n",
    "\n",
    "A system of nonlinear equations has the form\n",
    "\n",
    "$$\n",
    "  \\mathbf{f}(\\mathbf{x})=\\mathbf{0}.\n",
    "$$\n",
    "\n",
    "For a system, the derivative object is the **Jacobian matrix**. Each row of the Jacobian is the gradient of one residual component. That connection is the main reason the gradient material from Lesson 19 helps here.\n",
    "\n",
    "### Objectives\n",
    "\n",
    "By the end of this lesson, you should be able to\n",
    "\n",
    "- Use plots and signs to isolate roots of scalar nonlinear equations.\n",
    "- Implement and interpret bisection as a bracketing method.\n",
    "- Use Newton's method when derivative information is available.\n",
    "- Use finite-difference or complex-step derivative estimates when an analytic derivative is inconvenient.\n",
    "- Interpret a Jacobian as one gradient per residual equation in a nonlinear system.\n",
    "- Use `scipy.optimize.fsolve` for scalar equations and small nonlinear systems.\n",
    "\n",
    "### Key Terms\n",
    "\n",
    "- nonlinear equation\n",
    "- residual\n",
    "- root\n",
    "- bracket\n",
    "- bisection\n",
    "- Newton's method\n",
    "- secant method\n",
    "- convergence\n",
    "- derivative estimate\n",
    "- Jacobian matrix\n",
    "- nonlinear system\n",
    "- `scipy.optimize.fsolve`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear Equations\n",
    "\n",
    "Linear equations are characterized by *linear combinations* or the unknowns.  For example, the system of equations for $x$ and $y$\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    " ax + by &= 1 \\\\\n",
    " cx + dy &= 2 \n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "is linear because $x$ and $y$ appear with only constant coefficients.  There are no $x^2$ terms, or $\\sin(y)$ terms.  Just multiples of $x$ and $y$.  Any deviation from this pattern results in **nonlinear equations**.  For example, quadratic equation $ax^2 + bx + c = 0$ is nonlinear, and we know how to find the **roots** of this equations.  A root is the solution $x$ to any equation (generally, nonlinear) $f(x) = 0$.  Our goal now is to solve such equations using the tools at our disposal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytic Solutions\n",
    "\n",
    "In rare cases, nonlinear equations have closed-form solutions (i.e., you can write it down explicitly).  The quadratic equation is one example.  In fact, solutions exist for the general cubic (add a $x^3$ term) and quartic (add both $x^3$ and $x^4$).  The roots are hideous (compared to the quadratic equation roots), but they are known.  \n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**: Use SymPy to find the roots of the cubic equation $ax^3 + bx^2 + cx + d = 0$.  Can it also find roots to *quintic* equation $ax^5 + bx^4 + cx^3 + dx^2 + ex + f = 0$?\n",
    "\n",
    "*Solution*:  The cubic is easy, but SymPy can't do the quintic.  In fact, the roots to general polynomials of degree 5 and higher cannot be determined explicity, a fact [proven by Abel](https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem).\n",
    "\n",
    "<hr />\n",
    "\n",
    "For other problems, we often have to have a bit of luck.\n",
    "\n",
    "<hr />\n",
    "**Exercise**:  Find all values of $x$ such that $\\sin(ax) - b = 0$.\n",
    "\n",
    "*Solution*: The equation has no solution if $|b| > 1$ (because $-1 \\leq  \\sin(ax) \\leq 1$).  If $|b| \\leq 1$, then $x = \\frac{\\sin^{-1}b}{a} + n\\pi$ for $n = 0, \\pm 1, \\pm 2, \\ldots$.  Such an equation is called **transcendental** because it involves a transcendental function (i.e., $\\sin(x)$).  Transcendental functions cannot be defined in terms of polynomial of finite degree; in other words, if one needs an infinitely long [Taylor series](L20_NumericalIntegration.ipynb) to represent the function, it is transcendental.  Transcendental equations almost always require some form of numerical evaluation.\n",
    "\n",
    "**Exercise**: Solve the nonlinear system of equations $y - x = 4$ and $x^2 + y = 3$.\n",
    "\n",
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graphical Solutions\n",
    "\n",
    "When analytic methods fail, solutions  (albeit, approximate)  to nonlinear equations may be found **graphically**.  For single equation in $x$, it suffices to plot $f(x)$ over the range of interest and identify the roots.\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**:  Revisit $\\sin(ax) - b$ for $a=1$ and $b=1/4$ and find all roots between $x = -5$ and $x = 5$.\n",
    "\n",
    "*Solution*:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T18:55:46.423656Z",
     "iopub.status.busy": "2026-07-05T18:55:46.423572Z",
     "iopub.status.idle": "2026-07-05T18:55:46.705413Z",
     "shell.execute_reply": "2026-07-05T18:55:46.704908Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAEmCAYAAABrmPj4AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVE5JREFUeJzt3Qd0VEXbB/B/eiM9pJGQhJICoffekWJBEPFFURDxBURErIgKNhRfBD8LoIKiIIpKU0AgKr0mEDqEQBKSAGmk97bfmVkSAyQhQJK7e/f/O+cebpbZ7My9uzfPzp15xkij0WhARERERGQAjJWuABERERFRfWHwS0REREQGg8EvERERERkMBr9EREREZDAY/BIRERGRwWDwS0REREQGg8EvERERERkMBr9EREREZDBMla6AristLcWVK1dga2sLIyMjpatDRERERDcRa7ZlZWXB09MTxsbV9+0y+L0NEfh6e3vfrhgRERERKSwuLg5eXl7VlmHwexuix7fsYNrZ2aGuFRUVYfv27Rg8eDDMzMygRmpvI9un/3gO9Zvaz58htJHt039F9fwezczMlJ2VZXFbdRj83kbZUAcR+NZX8GttbS1fS40XNENoI9un/3gO9Zvaz58htJHt039FCr1HazJElRPeiIiIiMhgMPglIiIiIoPB4JeIiIiIDAaDXyIiIiIyGAx+iYiIiMhgMPglIiIiIoPBVGd6qqRUg7NXM3E0Ng3XsguRkVeEzPwiZOYVyf2iEg087C3h5WiFRg5W8HK0RiNHK/i52MDSzETp6hMR6b2iklKcupyBK+n5SMjMR1Km9t/EzHxkFxTD0doczjbmcG5gAecG2v2mDRugjbcDzEzY90SkFAa/eqK0VINzCVk4EHUNBy5ew+Hoa8jML672Ocfibn3M0swYPZs1xOAWbhgQ5CovykREVDMiwN0ZkYwdEUnYG5mCrILqr8OVaWBhiq5NnNGruYvcRKdETXKTElHtYPCr465lF2D1oVisOnQJiZkFt1xAO/k6yh5dO0sz2FtpNzsrMxgbGeFqRh4up+UhPi0Pl9PzEJeWi/TcIvx1NlFu4lrb0ccRg1q4YUTbRnC1s1SsnUREuio1pxA/HryEbWcScOpy5g3/52htJntz3ewt4WZrCXd7C7jZWcrrc1pukbyGi+enZIutACfi0+XjZddhQdyde7hdIzzZ3QeutrwOE9U1Br86KiIhCysPxWP9scsoLC6Vj9mYm6CTn5PsMejWxBktPe1gege3zjQabe/x9tOJCDmrvYiHxqTJbcH28xjbuTH+26cJPOyt6rBlRET6QXQgfLM7Gj8djkVeUUn542287NE3wBX9Al3RupE9jI2N7ugu3ukrmdhzIRl7zqfgyKU02TnxxY4L+HpPFEa2a4RnejVBM9cGddQqImLwq0PERfGfiGR8ecYY5w8cKH+8tZc9nu7hh2GtPGBuevfjxMRttSAPO7m9MLA5rqTnyZ6H9eGXER6bjhX7Y2Qv86OdvDClbzPZG0FEZGiiU3KwdOdFrAuPl/MnhOBGdniymy/6Bbiioe3dDxcTgXIrL3u5Te3bDLmFxXIYxTd7ouR1+OfQOLkNDHKTnRGdfJ1qsWVEJDD41SGFJaV4fd0ppOUaQ3QkDAl2l0FvBx/HOhkP5ulgJS/m47r6yHHEn/4dicPRqVh1MBZrQuPwSAdvvDzYn+OCicggpOcW4oPNZ7H2aDxKtTEvuvg54bl+zeTY3Lq4Dlubm8qODbGFxaTiq91R5UMixDa8lQfmPNiCwyGIahGDXx0isjBM6uWLsBMReHtsH/g2tKuX1xUX9O7NXOR2MOoaPvs7EvsvXpO3+rafTsD7I4IxtJVHvdSFiEgJW09dxZsbTstxucKAQFdM7dcUHXzqr+e1o6+T3C4mZ2PZnij8EhaPzSevYk9kMt4YFoQxnbw5MY6oFjD41TGTevqhUeZZxYYciPHEYguNScVbG07JMcJTfjyKB9p44t0HW8LRxlyRehER1YXkrALM+f0UtpxMkD+LsbbzR7WWd9yUIibQfTiyNR7v4oNZ607i5OUMvL7upByiNm9kK/n/RHT3mGiQKiXGmf0+rSee798MJsZG+OP4FQxatBvbTmv/QBAR6TMxAXh9eDwGLdolA19xnXuuX1Nser6nooFvRcGN7LF+ane8OTwIVmYmOBSdiqH/twdf7rgg54gQ0d1h8EtVEpPrXhocIC++zV0byNuB/115BDPXHENe4b8zn4mI9El+UQlm/nIcL645LtM/tvCww8bneuCV+wJ1bhEgkdFHZH/Y/mJv9PZvKLP//G9bBKasPobbpHonoiow+KXbau3lgD+e74kpfZvKiXjrwi9jzNcHZLJ3IiJ9Iq5bY74+KIcQiN7elwb5Y+O0HrKXVZd5O1nj+wmd8PEjrWXHhMgMtPCUicxMQUR3hsEv1YjoDXltSCB+mtRVJnU/EZ+Bh77ch9NXMngEiUgvnIzPwINf7MPxuHS5INAPT3fG8wOa681Sw2Jy8qMdvfHrf7vBzc4CiXlGGPXVIeyMSFK6akR6RT8+8aQzujRxxobneqBpQxtczcjH6KUHEHJGu0oREZGuEvMWRn+1HwmZ+fL6JYY59GjmAn3UxtsB6yd3hZ+tBln5xZiwIhRLdl6U45iJ6PYY/NId83G2wbqpPWTey9zCEjy7Mgxf7+aFl4h0jwgIF26PwPM/hSO/qBR9Axpi/XM94OtiA30mFtqY1qIEYzo2goh55289h5d+PY7iEu2KoERUNQa/dFfELcNvx3fCE10bywvvvC3n8NbGU5yBTEQ6Ffi+vfE0Pvvngvz52d5NsPypTrCzNIMaiAU/33uwBd4bESzHL687ehkvrDmGIgbARNVi8Et3TYyTe++hYMx9oIWcCCdWhnuTATAR6QCRCWzOH2ex8uAliIXZPhzZSi4UIYJENRHjgMUqnV+ObQ8zEyNsPnEV01YflVkhiKhyDH7pni+843v44ZNH28gAePWhWLz9+ymOPSMixYgcuL9GGeOn0HgZ+P7vkTb4T+fGqj4jQ4LdsfSJDjA3Mca204mYvOqITOlGRLdi8Eu14uF2Xlgwuo38QyN6gMWtRk6+ICIlAt+3fj+D/UnG8nr0yeg2eKSDl0GciAFBblj2VEdYiFRo55Iw6YcwBsBElWDwS7VmZHsv2cMi/uCIW41zfmcATET1G/i+tvYEfjlyGUbQ4H+jWsnrkiERC2F8N76TXBFuT2QKJnwXitxCroZBVBGDX6pVoofl41GtZQD8w4FLmMsAmIjqgbjT9Pq6E/j1SLwcgjWueSkeauNhkMe+ezMXfP90Z9iYm+BA1DU89+NRZoEgqoDBL9W60R29MX9ka7n//YFLWLA9gkeZiOrUopDz+CUsXk5oWzi6NTq4GHbO285+TvhhYmdYmhljR0QyZq/nXAyiMgx+qU482slbzq4WvtxxEb+ExvFIE1GdENeXsnRmH4wIxvBW7jzSADr4OOHz/7SXPeFrwuLw2d/aY0Rk6Bj8Up0Rs6uf799M7r+x/iT2RqbwaBNRrdp1Phmz1p+U+9P6NcNjKs/qcKcGtXDDuw8Fy/1Ff51nRwSRPga/ixcvhp+fHywtLdGhQwfs2bOnyrI7d+6Uqbhu3s6dO1evdTZkMwf546G2nigu1WDKqiM4n5ildJWISCVOX8nA1FVHUFKqwcPtGuGlwf5KV0knPdHVB8/1ayr3xReFHRFJSleJSFF6FfyuWbMGM2bMwOzZsxEeHo5evXph6NChiI2NrfZ5ERERuHr1avnWvHnzequzoRNfNj5+pDU6+zohq6BYzjxOzipQulpEpOeupOfh6RWhyCksQbcmzpgvJ9qqawGL2vTy4ACMbN9IflEQE+BOxmcoXSUixehV8Ltw4UJMnDgRzzzzDIKCgvDpp5/C29sbS5YsqfZ5rq6ucHd3L99MTEzqrc4EWJia4KtxHeDnYoPL6Xn474/hKGDudSK6S5n5RTLwTcwsQHPXBlg6rgPMxVq/VCXxxeCjka3Rs5kLcgtLMGFFqPwCQWSITKEnCgsLceTIEbz++us3PD548GDs37+/2ue2a9cO+fn5aNGiBd58803069evyrIFBQVyK5OZmSn/LSoqkltdK3uN+nit+tTA3AjfPNEOo78+hJOXM7Ey1xjDCguhRmo9h4bSPkNooz63T+Tynb46HOcSsuBqa4Fl49rB2vTGtuhz+2rqbtoo+sU/G9MaY5eHyuM3eVUYVk/sLBfF0DVqP4dqb58Sbbyjz4JGT5bhunLlCho1aoR9+/ahe/fu5Y/PmzcP33//vRzacDPx2O7du+XYYBHQrly5EkuXLpVjgXv37l3p68ydOxfvvPPOLY+vXr0a1tbWtdwqwxOVCXx5xgTFGiMM8y7BfV568fYjIh2xLd4IW+JMYGakwQvBJfBuoHSN9M+1fGDBCRPklhihu2spxjQtVbpKRPcsNzcXY8eORUZGBuzs7NTR81vm5jFdInavapxXQECA3Mp069YNcXFxWLBgQZXB76xZszBz5swben7F0ArRw3y7g1lb31xCQkIwaNAgmJmZQY3cwmLxxsZz+DPOBKP7d0DPZs5QE7WfQ7W3zxDaqK/tEyuW/XnwqNx/b0QwRrVvpKr23Yl7baNvqxRMXHlULgN9f/dWGN2h8mOpFLWfQ7W3T4k2lt2prwm9CX5dXFzkWN2EhIQbHk9KSoKbm1uNf0/Xrl2xatWqKv/fwsJCbjcTJ64+36D1/Xr1aXTHxvjjwBkcSDLGS7+dxKbne8LTwQpqo+ZzaAjtM4Q26lP74tNyMfO3kxD3Kv/T2RuPdfFVVfvu1t22sX8LD8wc6I9PQs5j7qazCPZyQGsvB+gatZ9DtbevPtt4J6+hewN9qmBubi6HL4hvERWJnysOg7gdkSXCw8Mwl7zUJaP8StHS0xapOYV4bvVRFBbzthsRVS6/qARTfzyK9NwitGpkjzkPtOShqgXP9WuGgUGu8vo7ZdVReT0mMgR6E/wKYjjCsmXL8O233+Ls2bN48cUXZZqzyZMnlw9ZePLJJ8vLi2wQGzZsQGRkJE6fPi3/f+3atZg2bZqCrSDBzBj4/LE2sLM0RXhsOuZtOcsDQ0SVeuePMzgRnwEHazMseaI9LM2Ysac2GBsb4ZNH28LX2Vpm4pn+U7hMhUakdnoV/I4ZM0YGtO+++y7atm0rJ7Nt2bIFPj4+8v9FDt+KOX9FhoiXX34ZrVu3ljmB9+7di82bN2PkyJEKtoLKeDtaY9GYtnJ/xf4Y/HH8Cg8OEd3g17A4/HQ4FmJqx/891g5ejpx4XJvsrczw1biOsDIzwd4LKfhk+62Tx4nURq+CX2Hq1KmIiYmR2RtE6rOKE9dWrFghMzmUefXVV3HhwgXk5eUhNTVVrgY3bNgwhWpOlRkQ5IapfbUrD72+9gQuJGXzQBGRdC4hE29uOCX3ZwzwRx//hjwydSDA3RbzH2kt9xfvvIh9F7gUPamb3gW/pM4lkMUKTWKlJrHykBjfR0SGTVwHXvjpGAqKS9E3oCGe799M6Sqp2oNtPOVEQmHmL8eQxvG/pGIMfklxpibG+Ow/7eDSwAIRiVn43zbediMydPO3npPXA5cG5lgwuo0cn0p16637W6CJi41cOe/1dSdkKlEiNWLwSzqhoa0FPn6kldxfvjcaeyN5243IUO0+n4zv9sXI/f890kZ+Maa6Z21uKsdVm5kYYdvpRKwJjeNhJ1Vi8Es6o3+gGx7v0ljuv/zrcaTnMu0OkaER6bZe+vW43B/X1Qf9Al2VrpJBaeVlj5cGB5Rn2YhK5jwMUh8Gv6RTZg8PkrfdEjLzMXvDKd52IzIg4ja7mPianFWApg1t8MawIKWrZJCe7dVEzsPIE+Oufz7GPOykOgx+Seduu4n0Z6bGRth84io2HLusdJWIqJ78EhaH7WcS5W13cfvdypz5fJUgxlcvHNNGpkE7eTkDi/46r0g9iOoKg1/SOW28HTB9QHO5//aG03JZUyJSt+iUHHmbXRC33YMb2StdJYPmYW+Fj0Zq52Es3XURBy5eU7pKRLWGwS/pJJH7t11jB2QVFOOlX45z1SEiFSsuKcWLa44ht7AEXZs4YVKvJkpXiQAMbeWBRzt6QSR9eHXtceQUFPO4kCow+CWdTX/26Zi2sDY3waHoVCzfG6V0lYiojogML8fi0mFraYqFj7aFCdOa6VT6M097S8Sl5jENJakGg1/SWT7ONvLCK3yy/TxiUnKUrhIR1TKRTWBhiHZM6VvDW8DTwYrHWIfYWprhw1Ha1d++PxCD0JhUpatEdM8Y/JJOe6yTN3o0c5arPL229gRKS5l0nUgtxOf59bUn5ee7V3MXjO7opXSVqBJiWenRHbTDH1777QRX4SS9x+CXdJqRkRE+fLg1rMy0wx9+ZtJ1ItVYdegSDsekyuFN8x5uJT/vpJvevL8FXG0tEJWSg0XXe+qJ9BWDX9J5jZ2t8fJ92qTr87acxdWMPKWrRET3SGRxmf/nObn/2pBAeDtZ85jqMJH2THxBEb7ZE4Xw2DSlq0R01xj8kl4Y390Xbb0dkF1QjNnrufgFkb4vZjFr3UnkFJagk6+jXMmNdN/AFm4Y0dYTYvTZq7+dQEFxidJVIrorDH5JL4jZ3x8/0lomv//nXBJ+P35F6SoR0V367Ug89kSmwNzUGB+Nai0XVSD9MOeBlnBpYI7IpGx8/vcFpatDdFcY/JLe8HezxfP9tYtfiGT417ILlK4SEd2hpMx8vLdJu5jFzEH+aNqwAY+hHnG0Mcd7DwXL/SW7LuLMlUylq0R0xxj8kl6Z3KcpAt1tkZpTiHev/wElIv0x5/fTyMwvRqtG9nimp5/S1aG7XPxiaLC7XHxo9oaTzMJDeofBL+kVcZt0vrhNagRsPHYFOyKSlK4SEdXQjnNJ+PNUghzGJD7HYjEb0t/hDzbmJgiPTWcWHtI7vPKQ3mnj7YAJPbQ9RnM2nmbOSSI9kFdYgrd/PyX3n+7hixaedkpXie6Bu70lXhqszcLz0Z9nkcJhaKRHGPySXnpxkD/c7CwQm5qLxTsvKl0dIrqNL3dckEvkethbYsZAfx4vFXiymw9aetrJYSzzNp9VujpENcbgl/RSAwtTvH1/S7m/dOdFRHPpYyKddSEpG1/t1n5JnfNAC9hYmCpdJaoFYtjKB3JxEmBd+GXsv5jC40p6gcEv6a1hrdzR278hCktK8fZG5v4l0tWcvm9tOIWiEg36BTTEfS3dla4S1SKRf/2JLto8zW9uOMXcv6QXGPyS3hJLob77YEs5CU7kDN188qrSVSKim4iJqQeirsHC1BjvPBjMJYxVSKzA6dLAAlHJOfh6V5TS1SG6LQa/pNd8XWwwpU9Tuf/uH2eQlV+kdJWI6LqMvCK8v1mbkvD5/s3kUuWkzqWP37o/SO5/vuMCYjgMjXQcg1/Se1P6NoWPszWSsgrw6V+RSleHiK5bsC0CKdmFaNLQBpN6N+FxUbEH23iiZzMXFBaX4u3fT8vhLkS6isEv6T1LMxO886B28tuK/TFccYhIB5yIT8eqQ5fk/vsPBcPC1ETpKlEdD0N7b0QwzE2Msft8MkLOJPJ4k85i8Euq0DfAVU6AEysOvcXJb0SKKi3VyJXcROffiLae6N7MhWfEAPi5iB5+bQ729zafYQ520lkMfkk13rq/BazNTXDkUho2HLusdHWIDNb68Mty5S/xeZw1TDsWlAzD1L7N4G5nKXM6L9vDyW+kmxj8kmp42FvhuX7N5P5Hf55DTkGx0lUiMjjZBcX4aOs5uf98/+Zws7NUukpUj0QO51nDAuX+lzsu4kp6Ho8/6RwGv6QqE3v6obGTNRIzC7B45wWlq0NkcD7/JxLJWQXwdbbG0z19la4OKTT5rbOvE/KKSjBvC1d+I93D4JdUN/lt9nDtbdZv9kQj9lqu0lUiMhhipcVv90aXD0PiJDfDnfw258EWMDYCNp24ioNR15SuEtENGPyS6gxu4VaecqcsxygR1b33Np2RK7n1DWiI/oGuPOQGrKWnPcZ2aSz35/5+GsUlpUpXiagcg19SZ6/DAy1gYmyE7WcSsTeS680T1bUd55Lwz7kkmBobyV5f8Tkkw/bSoAC5AMa5hCz8dDhW6eoQlWPwS6rU3M0W47pq15t/5w/2OhDVJXGXRfT6Ck/39EPThg14wAmONuZ4ebC/PBILtp9HWk4hjwrpBAa/pFovDvSHo7UZIpOyseqgNtk+EdW+FfujEZWSA5cGFnIZY6Iy/+ncGIHutnKp609CInhgSCcw+CXVsrc2w0uDA+T+wpDzSGWvA1GtE5kdPvtbm1nltSEBsLU041GmcqYmxph7fQXO1YdicT4xi0eHFMfgl1Tf6xDkYYfM/GIsZK8DUa0TXyxFbt/WXvYY1d6LR5hu0bWJM+5r6YZSDfDBZqY+I+Ux+CVVE5PexOQ34afDcYhkrwNRrTmXkIk1odqJTGKSm7HIbUVUiVlDg2BmYoRd55OxMyKJx4gUxeCXDKLXQaQ/KynVMOE6US3RaDSyF0/05g1r5Y5Ovk48tlQlXxcbPNVNu+iJeN8w9RkpicEvGYRZw4JkCqYdEcnYE5msdHWI9N7O8+KzlAJzE2O8NkS7nC1RdZ4f0Lx8EvLPoXE8WKQYBr9kEPxcbPBkhV4H0QtMRHdH9NqVjd0c38MXPs42PJR0WyLn74yB2tRni0LOIzO/iEeNFMHglwzG9AHNyhOu/xrGXgeiu/VTaBwuJGXLXrzn+jG1GdWcWPWtaUMbXMspxJc7tFlCiOqb3gW/ixcvhp+fHywtLdGhQwfs2bOn2vK7du2S5UT5Jk2aYOnSpfVWV9ItDtbmmD6geXnCdTFDnYjujOitE712wouD/OUXSqKaMjMxxuzhQXL/u70xiL2Wy4NH9U6vgt81a9ZgxowZmD17NsLDw9GrVy8MHToUsbGVL5sYHR2NYcOGyXKi/BtvvIHp06dj7dq19V530g1i1TdfZ2ukZBfgq10Xla4Okd4RvXUiZ7bovROpBInuVL8AV/Rq7oLCklLM33qOB5DqnSn0yMKFCzFx4kQ888wz8udPP/0U27Ztw5IlS/Dhhx/eUl708jZu3FiWE4KCghAWFoYFCxZg1KhR0FUm+flATg5gVkmPiokJYGn578+iXFWMjQErq7srm5srpnNXXtbICLC2vruyeXlAQUHVbbSxubFsaWnVda5YVvy+kpLbljU3NcYb/f3wwo9HsPKvMxjb0hkeDhXaLYj6inoLBQVAcTU9xJWVLSqqvH3i+IrjLBQWynJVupOy4v0g3hd3WlaUE+WrYmEBmJreWray9lUsK46BOBZVMTf/93l3UlacX/G6VRHlRPk7LSveY+K9VlHFNopzXF3ZisQxEMdCEJ8J8dmojbJ38rmvSdmy9om2VHyP3uYaEZenkb11wlv9fWCWn1c314iafu6rKlvWvopqeI2oUdl7vUbUxue+7POmi9eI25Q1KinBm/0aY+SZePxzJBpH2zVEe58K2ULKPm+COF7V/V6lrhF3W/bmz31Vf+uVvkbUVhyhqzR6oqCgQGNiYqJZt27dDY9Pnz5d07t370qf06tXL/n/FYnnm5qaagoLCyt9Tn5+viYjI6N8i4uLE1dsTUpKinxOXW85OTniLV7lVjJ06A3lS62tqy7bu/eNZV1cqi7bocONZX18qixbGhR0Y9mgoKrL+vjcUFa8TpVlXVxuLCvOa1Vlra1vLDt0aLXH7YayI0dWXzYt7d+y48ZVX/by5fKyxZMnV1/2/Pl/y86cWX3Z8PB/y775ZrVli/bv/7fshx9WXzYk5N+y//d/1ZfdsKG8bNGyZdWXXb3637KrV1dfdtmyf8tu2FBtWVHH8rIhIdWX/fDDf8vu31992Tff/Pc9ER5efdmZM/8te/589WUnT/637OXL1ZYV763ysmlp1ZcdOfKG93BtXSOKe/W6o2vElJWhGp/XNmnGfn1A568RRRYW8nqq1mtEbmioZsOGDbKNarxGiHaJ9uWtXavKa4Ro35bvv9f5a0TJPcQRJY0bl79H7zYmupNNxGkiXhOx2+3oTc9vSkoKSkpK4ObmdsPj4ueEhIRKnyMer6x8cXGx/H0eHh63PEf0IL/zzju3PL59+3ZYV+ydqEMPVfN/SUlJOLRlS/nPw0tKquy+T712DfsqlB1SWIjr3yFvkZGRgd0Vyg7KzUVVrc3KzsaOCmX7ZWfDroqyebm5CKlQtndGBhyrKFtYWIitFcr2uHYNLlWUFe+FLRXKdklKgjuqVrFsx4QENKqmrLibUHL9W3G7+HhUd2P3r7/+QqG9vdxvfekS/Kopu2PHDuRdfz+2iIqCdvRx5cRY9qxLl+R+QGQkqksktW/fPqQnaZPGNzt3DtqFRCt38OBBXLv+zd3v9Gm0rqasuEuSeH3f+/hxtK+mrBhWdOX658MzPBydqil74vhxxF0/H25hYehaTdnTp08j+npZ55Mn0bOasufOncOF62UdIiPRp5qykZGRiLhe1jY2Fv2rKRsVFYUz18taJSZicDVlYy9dwonrZc0zMjC0mrLx8fEIv15W9FLeX03ZqwkJCKvwHq6ta0RaamqNrxHJqRnYcioRRtCgh02i/Gzr8jVCCAkJUe014sCBA0DjxrKNqrxGuGjP7LHwcPVeI1A9XbhGpN5LHHG9B7zi57Au5d5BT7ORiIChB65cuYJGjRph//796NatW/njH3zwAVauXCnf1Dfz9/fHhAkTMGvWrBsuAj179sTVq1fh7n7rpbCgoEBuZTIzM+Ht7S2DZTu7qi7ftaeoqAj//PEH+vfvDzOVDnsoKijAP//8U3kb63jYQ8Wyb6w/jU0nE9DJ1wHLxrWHUdmtyXu8pSnPYWXtU8mwh0rbp7JhDze0UYXDHsrbN3AgzCpe16r4veLPxPgfwrH/ah5GtffERw8H1+3QqHsc9lDevgce+Pc9qrJhD0Wmpgj5+28MGjQIZuLY6tA14rZlK3zuF++6iKW7YuDtZIn1k7vJoWni81Z0PWga1K8fzKp7P+jpsAfxHg3Zvh2DevSo/G99hbJKXCNqI44oyshAyL592vdoVW2sRSJec3FxkUH47eI1ven5FQ0yMTG5pZdXfIO5uXe3jAhuKytvamoKZ2fnSp9jYWEht5uJE1cfJ08QPQpmDg41ez0Hh5r/4jspe72notbLijYVFdWsjXdyvO+i7Asj2mHjxV3YfSUf+xKL0C/QtXbqUJP21XHbaly2pnczKpa9XfvEYxUvlrf7vXdStuIFu7bKCjd/5qtrYyXXhypVHLtYm2Xv9XNf1j47uxvbV8Xv3X46QQa+FqbGePm+QO1z6vIaca9ly9pX8bqtK5+52ip7PYC9479N9XGNuIPP/TND2mDlyTSczy7ArxEZGN/D78b2WVnVvH31eY2ojbJGRjX/W1/f14hajCPqK366k9fQm2wP5ubmMmXZzd3n4ufu3btX+hzRQ3xzeTF8oWPHjvUWyJLu8nK0xoTu2oUvPvyTC18QVbWgxUfXZ+RP7OkHD/saflkhqgEbC1O8OEg7wOP//o7kwhdUL/Qm+BVmzpyJZcuW4dtvv8XZs2fx4osvyjRnkydPlv8vhjc8+eST5eXF45cuXZLPE+XF85YvX46XX35ZwVaQLpnaV7vwxfnEbKw9Eq90dYh0zpqwOEQl58DJxhyT+zZVujqkQmM6esvUeWm5RViykykoqe7pVfA7ZswYmbbs3XffRdu2bbF79245ScHHx0f+vxjHWzHnr1gMQ/z/zp07Zfn33nsPn332mU6nOaP6ZW9thuf7a1eo+iQkAnmF1YzzIzIwYiGYRSGRcn96/2aws+QdM6p9pibGeH2oduGLb/dG40p6NWNmiWqB3oz5LTN16lS5VWbFihW3PNanTx8cPXq0HmpG+mpcNx+s2B+D+LQ8fLsvmsu1El33ze4ouSCMj7M1xnbRdjIQ1YWBQa7o7OeEw9Gp+GT7eXz0cAseaKozetXzS1QXLExN8Mp9AXJf3HK7ll1N9gEiA5GUmY9v9kTJ/VfvC9TOwieqIyLbzuxh2t7fdeHxOHs1i8ea6gyvZkQAHmjtiVaN7OVt3s//ucBjQgZv0V+RyC0sQVtvBwxrVV2GXKLa0cbbAQ+08ZRZveZvO8/DSnWGwS+RTE9ohFnDtGniVx28hJiUanIZEqnchaQsrAnVzp94Y1jQvzmwierYK4MDYGZihH0Xr+FcOt93VDcY/BJd172pC/oFNERxqQYfb7t10RQiQzF/awRKNcCgFm5yHCZRfWnsbI1xXbUpKP+INUapeCMS1TIGv0QViBnHxkbAlpMJCI9N47EhgxMWk4qQM4nyc/DaEO1YeKL6NK1/M9hYmCA+xwibT924UBVRbWDwS1RBgLstRrX3kvsf/XlOLutKZCjE+12874UxnbzRzNVW6SqRARI5pZ/tqV3pbeFfF1BYXM3yxkR3gcEv0U1eHOQvZ7Yfik7FzvPJPD5kMESPb9ilNFiaGeOFAf5KV4cM2PjujWFnppEpKFcfuqR0dUhlGPwS3cTTwap82eP5f55DCceckYEsY/zxtgi5/3QPP7jbWypdJTJg1uamGOqt7fH97J8LyMovUrpKpCIMfokqMaVvU9hZmuJcQhY2HrvMY0Sqt/ZoPC4kZcPB2gz/7cNljEl5XVw1aOJijdScQrngClFtYfBLVAkHa3NM7Xd92ePt55FfxGWPSb3Est5lyxhP69cM9lZcxpiUZ2IEvDSoudz/Zk80krLyla4SqQSDX6IqjO/uC3c7S1xOz5O5f4nUSizvnZCZj0YOVniiK5cxJt0xKMgV7Ro7IK+oBP/3l/YLGlG9B78RERGYO3cuBgwYgKZNm8LDwwOtW7fGU089hdWrV6OggEvDkjpYmplg5iDtpJ8vdlxAJseckQql5xZi8U7tqoYvDfaX73siXSEWWJk1VLvs8c+hcYhKzla6SmRIwW94eDgGDRqENm3aYPfu3ejUqRNmzJiB9957D0888YRMkTN79mx4enpi/vz5DIJJFUa2b4Tmrg2QnluEr3ZdVLo6RLXuyx1iMlExAt1t8VDbRjzCpHPEQisDg1zl5OP/XZ+USXQvTGtacMSIEXjllVewZs0aODlVveLPgQMHsGjRInzyySd444037qlyREozNTHGq0MCMemHMCzfG40nu/nCzY6z4EkdrqTn4fsD2iE9rw0NhIlY2YJIB4nr8D/nkvDnqQQci0tHW28HpatEhhD8RkZGwtzc/LblunXrJrfCwsJ7rRuRThA9Dp18HREak4ZP/zqPD0e2VrpKRLXi//65KBcQ6NrECX39G/Koks7yd7PFyPZe+O1IPD768yx+mtRVDokgqtNhDzUJfIXc3Nw7Kk+k68QF9vWhgXL/lzBtOigifXclF1h/7Er5st4MJEhfFiA6GJWKXVyAiOo720Pfvn0RHx9/y+OHDh1C27Zt76U+RDqpg48TBrVwk2POFnDMGanA5lhjiNW7hwa78xYy6QWRjeSpbtpsJPO3RqCUCxBRfQa/dnZ2MsPDzz//LH8uLS2VGSB69+6NBx988G7rQqTTXr0vAGJI5NbTCTgam6Z0dYjumljC+FSasRzj+/J9ATySpDem9m0GWwtTnL2aid+Pa+9cENVL8Pv7779j3rx5eOaZZzB27Fj07NkTy5Ytw+bNm7FgwYK7+ZVEOq+5my0e6eAl9z/685zMcEKkb8T79n/btflSH2nfCE0bNlC6SkQ15mhjjsl9tSsQLtgegYJiLkBE9bjIxeTJk/H888/L3t+wsDD88ssvGDhw4N3+OiK9MGOgdszZ4ehU7IxIVro6RHcs5Ewijsamw8xYg+f7NeERJL3zdA8/uNpaID4tD6sPxSpdHTKU4DctLQ2jRo3CkiVL8NVXX+HRRx/F4MGDsXjx4tqvIZEO8XSwkiu/CfO3npNjgIn0RcU8qX08NEzbR3rJytxEdkQIn/8j8lQXKV0lMoTgNzg4GImJiXLhi0mTJmHVqlVYvnw53nrrLQwfPrz2a0mkQ6b2bQpbS1OcS8jCxmOXla4OUY2tPRqPyKRsOFiZYYBnKY8c6a1HO3qhiYsNUnMK8c2eaKWrQ4YQ/IohD2KVNz8/v/LHxowZg+PHjzO/L6meg7U5plwfc/bJ9vMcc0Z6Ib+oBItCzsv9yX38YF3jLO9EurkA0SvXJ2su2xOF5KwCpatEag9+RQ+vsfGtT/Xy8kJISEht1ItIp03o7gc3OwtcTs/DqoMcc0a674cDMbiakQ9Pe0s80dlb6eoQ3bMhwe5o4+2A3MISfP6PdhInUa0Gv7Gxd/YH/vJl3g4mwxhz9sU/kcjkmDPSYRm5Rfhyx8XyhQIszEyUrhJR7SxANES7AJGY+BaTksOjSrUb/Hbq1EmO7z18+HCVZTIyMvDNN9/IMcHr1q2r6a8m0kujO3ihSUMbpOUW4etdUUpXh6hKS3ZdREZeEfzdGsglYonUoltTZ/QNaIhisQDRdu1kTqLbqfGor7Nnz8rcvkOGDIGZmRk6duwIT09PWFpayuwPZ86cwenTp+Xj//vf/zB06NCa/moivR1z9up9gZi86giW743Gk9184GpnqXS1iG5wNSMP3+3TTgh6bUigXNiilKlRSUXEdVgsd7zpxFU82zsdrb0clK4SqaXnVyxnPH/+fFy5cgVLly6Fv78/UlJSEBmpHWfz+OOP48iRI9i3bx8DXzIY97V0Q7vGDsgrKsGnf3PMGemeT0MiUVBcis5+Tugf6Kp0dYhqXQtPOzzctpHc5wJEVKs9v+3atUNCQgIaNmyIl156CaGhoXB2dq7p04lUPeZszNcHsSY0DhN7+qGxg4XS1SKSIhOz8OuROLn/+tBA+X4lUiMxll30/O6/eA27I1PQx7+h0lUiNfT8Ojg4ICpKO64xJiYGpaXMEUkkdGnijAGBrtoFBLZyzBnpjvlbIyDWYRnS0h3tGzsqXR2iOuPtZI1x3XzKe39LuQAR1UbPr1jRrU+fPvDw8JC9B2Jsr4lJ5TOGy4JkIkPx6pBA7IhIwtbTCQiPa6x0dYgQGpOKv84myjG+rwzR5kMlUrPn+jXDL6FxOHs1E78fv4IR7bRDIYjuOvj9+uuvMXLkSFy4cAHTp0+XmR9sbW1r+nQiVQtwt8Wo9l749Ug8Pt52Hk94KF0jMmQajUb2fgmPdvRG04YNlK4SUZ1zsjHH5L5N5RLeIvPD0FbusDBlWj+61R2t8SMyPQhiYtsLL7zA4JfopjFnorch7FI62lkZgQt9k1JCziTiyKU0WJoZY8bA5jwRZDAm9PDF9/tjEJ+Whx8PxuLpnv+uREt0Tyu8fffddwx8iW7i6WCF8T185f4fl4zlGGCi+lZcUoqPt2nHnosJmG5Mv0cGxNrctHwBIrHqGxcgoloLfomoclP7NIOdpSkS8oyw4dgVHiaqd78diceFpGw4Wpvhv32a8gyQwXm0478LEH21S7uyIVFFDH6JapG9tRkm99HeZvv07wvIL+JqAlR/cguLsTDkfPnkHztLMx5+MsgFiMqWPV62J1ou9EJUEYNfolo2rktjOJhrkJBZgO/2xfD4Ur1ZvicaSVkF8HayKk/7RGSIBrVwQydfR7nAy6LrXwiJyjD4JapllmYmGN5Ymwd78Y4LSM0p5DGmOpecVYCl12/xvnJfIGe5k0ETKVlnDQsqHwp0LiFT6SqRDmHwS1QHOrpoEOhui6yCYjnpgqiuffZ3JHIKS9Dayx73t2KuPSKxsMuwVu5yoZf511P/EQkMfonqgLER8Np92hnHqw5ewqVrOTzOVGcuJmdj9eFYuT9raBCMxRuQiORdEFNjI+yISMb+Cyk8IiQx+CWqIz2bOaO3f0MUlWjwMZc9pjr08dZzMrWeWGa7W1NnHmui6/xcbPB4F+2qmx9y2WPSt+A3LS0N48aNg729vdzEfnp6erXPGT9+vBz3U3Hr2rVrvdWZaNbQQBgZAZtPXkV4bBoPCNW6sJhUbDudKO82vD5UO8OdiP71/IDmaGBhipOXM/DHCaagJD0KfseOHYtjx45h69atchP7IgCuyap0V69eLd+2bNlSL/UlEoI87OSyx8K8LWflsrNEtUW8n8T7ShjTyRvN3bjkPNHNXBpYYHKfJnJfLH1cUMwUlIZOL4Lfs2fPyoB32bJl6Natm9y++eYbbNq0CRER2pWMqmJhYQF3d/fyzcnJqd7qTSS8NNgfFqbGCI1Jk8vOEtWWracScDQ2HVZmJuWrWhHRrSb2bAI3Owu57PHKA5d4iAycKfTAgQMH5FCHLl26lD8mhi+Ix/bv34+AgIAqn7tz5064urrCwcEBffr0wQcffCB/rkpBQYHcymRmatOjFBUVya2ulb1GfbyWUtTexpvb52JtigndfbB0dzQ++vMsejZ1hJmJXnzvNMjzpy9tLCopxUfXZ7A/3cMHTlYmNa6vPrTvXqi9fYbQxtpun6kR8EL/pnhjwxmZgeeh1u5wsFZuERi1nz8l2ngnr2Ok0YP7sPPmzcOKFStw/vyNiar9/f0xYcIEzJo1q9LnrVmzBg0aNICPjw+io6Px1ltvobi4GEeOHJE9wpWZO3cu3nnnnVseX716NaytrWupRWRo8oqB98JNkFNshNF+JejprvMfO9Jxu64aYV2MCRqYafBWuxJYmihdIyLdJlKefXzCBFdzjdDHoxQjfbX52EkdcnNz5RDZjIwM2NnZ6W7Pb1WBZkWhoaHyXzFZ7WYibq/s8TJjxowp3w8ODkbHjh1lILx582aMHDmy0ueIQHrmzJk39Px6e3tj8ODBtz2YtfXNJSQkBIMGDYKZmTqXJlV7G6tqX65rLN7dfA5/J1li1thesLXUixsvBnf+9KGN6blFePvTPQCK8erQFhjZyVtV7btXam+fIbSxrtrnEHANE74/gv1JJnhzTC/4OttACWo/f0q0sexOfU0o+td32rRpeOyxx6ot4+vrixMnTiAx8daxksnJyXBzc6vx63l4eMjgNzKy6kUHRI9wZb3C4sTV5xu0vl9PCWpv483tG9fdD6sOxSEqJQff7LuE166vPa+v1H7+dLmNS/dEIiOvGP5uDTC2iy9M73IYja62r7aovX2G0Mbabl+/IHf0DWiInRHJWBByAV+N6wglqf381Wcb7+Q1FA1+XVxc5HY7YoKb6MY+fPgwOnfuLB87dOiQfKx79+41fr1r164hLi5OBsFE9U2M831jWBCe+SEMy/dGY2znxvB24lAaujPRKTn44UCM3J89vMVdB75Ehmr2sCDsiUyRKQIPRl1D1ybMjW1o9OKqGRQUJFOWTZo0CQcPHpSb2L///vtvmOwWGBiI9evXy/3s7Gy8/PLLcrJcTEyMnPj2wAMPyGD74YcfVrA1ZMgGBLmie1NnFBaXYv5WLrdJd05MmhQLp/Txbyg3IrozIiXgfzprhwp9sPksSsVgYDIoehH8Cj/++CNatWolx96KrXXr1li5cuUNZUTaM9EbLJiYmODkyZN46KGH5MS4p556Sv4rgmFbW+bCJGWIMepvDm8hF77YdOIqjlxK5amgGjsUdU32VpkYG2H28CAeOaK7JFIDli18seHYZR5HA6M3M25Eft5Vq1ZVW6Zi4gorKyts27atHmpGdGdaeNrh0Q7eWBMWh3c3ncX6Kd1hLJbnIqqG6J16f7N2QQvRa+XPBS2I7mnhi+f6NZN34MTCF0ODPWBlzpQphkJven6J1OSl+/xhbW6C43HpXG6TakT0ToleKlsLUy5oQVQLJvTwRSMHK1zNyMeyPVE8pgaEwS+RAlxtLTG1b1O5P//Pc8gv4nKbVLW8whJ8vFW7muVz/ZvJXisiujeWZiZ4bag2686SXReRlJnPQ2ogGPwSKeSZXk3gaW+JKxn5MvsDUVW+2ROFhMx8eDlaYXx3Xx4oolryQGsPtPV2QG5hiRz+QIaBwS+RDvQ6LN5xAUlZ7HWgW13NyMOSnRfl/utDA+X7hohqbxLy2w+0kPu/HomXQ9FI/Rj8EinowTaestchp8JtbaKKPvrzHPKKStDRxxHDWzFHOVFta9/YESPbNZL7c/84zdRnBoDBL5GO9Dr8diQex9jrQBWExqRi47ErMjXe3AdbVrucOxHdPXEXTkxCDo9NZ+ozA8Dgl0gHeh1GtfeS+3N+Z68DaZWUajD399Ny/7FOjRHcyJ6HhqiOuNlZYlr/ZuV3W3IKinmsVYzBL5EOeG1IgEy4LsabrT0ar3R1SAesCY3D6SuZsLU0xcuD/ZWuDpHqTezpBx9nayRlFeDLHReUrg7VIQa/RDrA1c4S0wdoex3mb41AZn6R0lUiBWXkFmHBdu0Y8JmD/OHM1GZEdc7C1ASzh2lXTly2JxqXruXwqKsUg18iHTG+ux+auNggJbsAn/0VqXR1SEGL/jqP1JxC+Ls1wBNdfXguiOrJoBZu6NXcBYUlpeUrKpL6MPgl0hHmpsblk99W7I/BhaQspatECjifmIWVBy/J/TkPtISZCS/TRPU6Cfn+FjAxNkLImUTsiUzmwVchXlWJdEjfAFcMDHJDcakG7/xxBhqNRukqUT0S5/udP07LyW5DWrqjRzMXHn+ietbczRZPdtPecRHX4aKSUp4DlWHwS6Rj3ro/COYmxtgTmSJ7HshwbDudgH0Xrsm7ALOHa8ceElH9mzHQH0425riQlI0V+2J4ClSGwS+RjvFxtsGk3n5y/73NZ5BfVKJ0lageiNRK7/5xRu7/t3cTeDtZ87gTKcTeykyuqFg2Bv9Keh7PhYow+CXSQVP7NoO7nSXiUvOw+PrStqRun/0diSsZ+fBytJLnn4iU9Uh7L7myYm5hCd7bpP1iSurA4JdIB9lYmGLO9clvS3dexMXkbKWrRHUoIiELy/dGy/13H2oJK3MTHm8ihRkbG+G9EcFy8tufpxKwIyJJ6SpRLWHwS6SjhgS7o19AQ5ly560Npzj5TaVKSzV4c8NJOcnxvpZu6B/opnSViOi6IA87TOjuK/fnbDzNYWgqweCXSIdT7rz7UDAsTI2x/+I1bDx2RekqUR347Wg8QmPSYG1uIlObEZFumTHIXw5Di03NxRIOQ1MFBr9EOkxMepo+oLncf3/zGbnyF6lHWk4hPtyiTaQ/Y2BzeDpYKV0lIrqJWHq+LAe7CH6jU7jym75j8Euk4yb1aoJmrg2Qkl2Ij7edU7o6VIvmbz2HtNwiBLjZYkIPbYYPItI9Q4Pd0dtfOwzt7Y0chqbvGPwS6TiR8/WDEcFyf/XhWByNTVO6SlQLwmJS8XNonNz/4OFgruRGpOvD0B5sKa/HIgf7lpMJSleJ7gGDXyI90KWJMx7p4AWx4Nvs9adQzBWH9JpYMerNDafk/qMdvdDR10npKhHRbfi62GBKn6ZyX6zEmJHHYWj6isEvkZ6YNTQQDtZmOHs1Eyv2c8UhfbZsTzTOJWTJ8/n6UK7kRqQvpvRtiiYuNkjKKigfr0/6h8EvkZ5wbmAhA2Dhk+3nEZeaq3SV6C6InM1ixSjhzeEt5BKqRKQfLM1M8NGo1nJfDFvadyFF6SrRXWDwS6RHRnfwRtcmTsgrKsFra08w968e5vR97bcTKCwuRR//hhjVvpHSVSKiO9TZzwnjuvrI/VnrTiK3sJjHUM8w+CXSsxWH5o9qDUszbe7fnw5rJ0yRfvjhQAzCLqXBxtwE80a2kpNoiEj/vDokAJ722ty/C7dr7+SQ/mDwS6RnfJxt8Mp92uEP87acxeX0PKWrRDUghql8vC1C7r8+LAiNmNOXSG/ZWprhg5Gt5P63+6IRziw8eoXBL5EeGt/dFx18HJFdUIw31p3k8Acdp9Fort8eLZG3TB/v3FjpKhHRPeoX4IqR7RqhVAM5DE0MZyL9wOCXSA+ZXB/+IHJO7jqfjN+OxCtdJarGL2Fx2HshRS5VLc6bGL5CRPrvrftbwNnGHOcTs7F45wWlq0M1xOCXSE+JVd9mDvKX++9tOoPEzHylq0SVEOfl/c3alEgvDfaHn4sNjxORSjjamOOdh1rK/S93XEBEQpbSVaIaYPBLpMee6emHNl72yMwvxuz1HP6gi8MdxKIkWfnFaOPtgIk9myhdJSKqZcNbeWBQCzcUlWjw4ppjHP6gBxj8EukxUxNjfPxIG5iZGOGvs0nYeOyK0lWim4Y7/HU2UZ6fj0e1lsNViEhdRNYWsQS9o7UZzlzNLM/jTbqLwS+Rngtwt8X0/s3l/lsbTyE+jYtf6ILolBy888cZuf/S4AB5nohInVztLPHhSO3iF0t3XcTh6FSlq0TVYPBLpJIlN9s3dpC318VttxIx/ZgUU1RSihlrjsnsDmJRkkm9ONyBSO2GBLtjdAcvaDSQ1+HM/CKlq0RVYPBLpJLhD5+OaYcGFqYIjUnD4h2cdaykz/+OxPG4dNhZmmLho2053IHIQMx5sCW8naxk/vW5v59WujpUBQa/RCrR2Nka716fdfzp35FMuq6QI5dS8cX1Lx8fPNwKnlzMgshgiA6IRY+2hRjev+7oZWw+cVXpKlElGPwSqcjD7RrhwTaectjDCz8fk4tgUP3Jyi+Swx3EqBOR/P6BNp48/EQGpqOvE6b2bSb331h/EgkZTEOpaxj8Eqls1vF7I4Ll0rlizXnedqtfc38/g7jUPHg5WpXn/iQiw/PCwOZo1cgeGXlFeOW34yjlPAydwuCXSGXsrcywaIz2tptY+W3TCaY/qw/i9ubao/HyuIvjb2tpVi+vS0S6x8zEWF4HLM2MsScyBUt3X1S6SlQBg18iFers54Tn+mlvu81adxJxqUx/VpeikrPx2toTcl/c7uzk61Snr0dE+rEK55wHtHeAFmyLwP6LKUpXia5j8EukUtMHNEdbb236syk/HkF+UYnSVVKl3MJiTFl1VI6v7uzrJG93EhEJj3Xyxsj2jeQ8gOk/hXMZeh3B4JdIxbfdvny8PZxszHHqcibe2nBKLrdLtUccT9GzHpGYhYa2FvhibDt53ImI/l39rRUC3W2Rkl2I51eHyzzgpCy9uUp/8MEH6N69O6ytreHg4FDjP0xz586Fp6cnrKys0LdvX5w+zbx7ZDjExLfP/9NOjkP99Ug8Vh+OVbpKqvLDgUtySWmxbPGXY9vLVZ6IiCqyMjfB4sfbyzRoh2NS8b9tETxACtOb4LewsBCjR4/GlClTavycjz/+GAsXLsQXX3yB0NBQuLu7Y9CgQcjKyqrTuhLpkh7NXPDKfYFyX2R/CI9NU7pKqnDkUhre36xdvnjW0EA5zpqIqDJNGjbAgtHa5Y+/3h2F7WcSeaAUpDfB7zvvvIMXX3wRrVq1qnGv76efforZs2dj5MiRCA4Oxvfff4/c3FysXr26zutLpEsm92mCIS3dUVSikeNTk7MKlK6SXkvJLsBzPx6Vx3N4aw9M7OmndJWISMcNCfbAM9evFa+tO42kPKVrZLhMoVLR0dFISEjA4MGDyx+zsLBAnz59sH//fvz3v/+t9HkFBQVyK5OZmSn/LSoqkltdK3uN+ngtpai9jbravnkjWuB8YhaiUnIwbfURrHiqg1wWWS3tq03VtbG4pBTTfjyChMx8NHGxwfsPBqG4WL8WE1H7OVR7+wyhjWpt38yBTXEsLg1hl9LxbYQJHsrOg2MDqFJRPZ/DO3kdI42ezYBZsWIFZsyYgfT09GrLiQC3R48euHz5shzzW+bZZ5/FpUuXsG3btkqfJ8YIi17mm4neYjHemEifJeQCC0+aoKDUCH09SvGwLyde3AlxtVwbbYw9icawMNZgZqsSuPOyQER3IKMQ+N8JE2QVGSHIoRSTAkthYsRDeK/Enf2xY8ciIyMDdnZ2utvzW1WgWZEYq9uxY8d7mmlZkYj1b36solmzZmHmzJk39Px6e3vLHuTbHcza+uYSEhIixyabmakzSb7a26jr7fMITMD0NSew86ox+nZsicc7e6uqfbWhqjYu3xeDPYnnIS4hnzzaFve1dIM+Uvs5VHv7DKGNam9fQNtreOK7MJxNN8bhksZ494GgamMTfVRUz+ew7E59TSga/E6bNg2PPfZYtWV8fX3v6neLyW2CGPrg4eFR/nhSUhLc3Kr+gyWGRojtZuLE1ecHsL5fTwlqb6Outu/Bdt6IvpaPRX+dx7ubzsLTwRqDW2o/L2poX22q2EaxgttHW8/L/dnDgnB/Wy/oO7WfQ7W3zxDaqNb2tfd1xlPNS7H8vAl+Do2XE+Ke7d0UamRWT+fwTl5D0eDXxcVFbnXBz89PBsDiW0e7du3KM0bs2rUL8+fPr5PXJNIX0wc0w9WMPPwcGofpP4dj9aSuaN/YUelq6aywmFS8+MsxuT++uy8nuBHRPWvlpMGsIQGY92cE5m05B29Hawxt9W9nHdUdvcn2EBsbi2PHjsl/S0pK5L7YsrOzy8sEBgZi/fr1cl/cPhBjg+fNmycfO3XqFMaPHy/H7YoxIUSGTHw+3h8RjH4BDZFfVIqJK0LlEr10K3FcnvkhDIXFpRjUwg1v3d9CdbcniUgZ47s1xlPdfOT+jDXHcJSpKOuF3gS/b7/9tuzBnTNnjgx4xb7YwsLCystERETIgc5lXn31VRkAT506VY4bFpPftm/fDltbW4VaQaQ7RKaHL8a2R2sve6TlFmH8d6FMgXaTa9kF8rik5xahjbcDPnusnVzQgoioNogv0uILdf9AVxQUl2LS92GIvZbLg1vHjPUpy4OYrHbzJlZtKyN+Fr27Fd9UYlLd1atXkZ+fL4c8iHy/RKRlY2GKb8d3QmMna8Sm5mLi96HIKdCvtF11Jb8YePbHcHlcxPFZ/lRHuVITEVFtd0SIlThbetrhWk4hHl9+EJfTmQS4LulN8EtEdcOlgQW+f7ozHK3NcCI+A5N+CENuoWEHwFn5RVhy1gQn4jPhYG2GFRM6yeNERFSXHRG+ztaIS83DY18fwBUGwHWGwS8Rwc/FRl54bcxNsP/iNYz/NhTZBtoDnJFXhPHfH0FMthHsrUyxamIXORObiKguudlZ4qdnu8KnPAA+yAC4jjD4JSKpXWNHrHymC2wtTHE4JhVPLj+EzHx1ra5Uk8BXtFv0+FqbavD9+I4IbmSvdLWIyEB42Fvhp0ldy4ei/eebgzIzD9UuBr9EVE6kO/txUhfYW5nhaGw6nlh2CBm5hhEAp+cWyvYej8+QQ0CmtSiRY/CIiOqTp4OV7AH2drLCpWu5+M/XB5GQkc+TUIsY/BLRDVp7OWD1pC5wsjGXY4BFz0NqTqHqA9/Hlx3CycsZst0/TOiIRjZK14qIDFUjByv8/Gw3GQDHiAD4m4OIT2MWiNrC4JeIbtHS017eehOTvM5czVR1z8OlazkYvfQATl/JhLONuWx3oDvTIRKR8gGwuB55OVohOiUHI77chyOX0nhaagGDXyKqVIC7LX5+titcbS0QkZiFB77YiyOXUlV1tPZfTMFDX+5DZFI23Ows5K1G0W4iIl3g5WiNXyd3QwsPO6RkF8oe4I3HLitdLb3H4JeIqtTMtQHWTukue0KTswrk7ONfwuJVccRWHbyEJ5cf1i5g4WWP36f1hL8bA18i0r1JcCIAHhjkJleafOHnY1gUcl6ubUB3h8EvEVXL28laBsBDg91RVKLB7I1n8GuUMYpKSvXyyIl6v7XhFN7ccArFpRo81NYTa/7bTaYZIiLS1TzAX43rgP/2biJ//r+/IzH952PILypRump6icEvEdXowrv48fZ4ebA/jIyAvYnGeGrFEaRkF+jdcsVPfXsYKw9ekj+/cl8APh3TFpZmXLmNiHSbWFp91rAgzB/VCqbGRvjj+BU8+tUBXEjKVrpqeofBLxHViFgufFr/5lgyti0sTDQIjUnD8M/2YPvpBJ0/guL2oBgnN3DhLrmIh7W5Cb4e1wHP9Wsm20VEpC/GdGqMlRO7yNUnRUaeYZ/twTe7o1BSymEQNcXgl4juyIBAV8wMLkETF2skZhbg2ZVHMGXVESRm6mY2CJGlQizZLMbJpeUWyfHL66Z2x+CW7kpXjYjornRr6ow/X+iFPv4N5TjgD7aclb3AUcnsBa4JBr9EdMfcrYGNU7thSt+m8lbcn6cSZK/qj4cuoVRHeh9Eb+/Ph2MxaNEu/HU2CWYmRpg5yF9ObAt05+IVRKT/E+FWTOgkh0E0sDCVadCG/t8eLNvDXuDbYfBLRHdFjJN9bUgg/pjWU2ZLyMovxuz1pzDm6wM4dTlD0aMaGpOKsd8cwuvrTsp6tfF2wObpvTB9QHOYm/KyR0TqIIZtiWEQ217sjZ7NXFBQXIr3N5+VQ9K2nLyqaGdEYmY+4nOgk0yVrgAR6bcWnnZYN7UHvt8fgwXbI+RY4Ps/3ysvxM/2boJezV3qZVyt6OndeT4ZS3ZcxOEYbT5iC1NjOaltQg8/2UNNRKTWBTFWTuyM1Ydj8eGWcziXkIWpPx6Fv1sDPN+/OYa18qi3a+CJ+HR8uzcam05chbuVCSbpYEo2Br9EdM/ERfXpnn64L9gd8/88h80nr2LvhRS5iTG2Igh+oI0nzExqv9dVTPL489RVLN5xUa5GJ5ibGGNUh0aY2reZTNVGRKR2opPh8S4+GN7KA9/ui8F3+6JxPjEbz/8ULlOjPdevKYYGe9RJdpviklJsP5Mog96wCqvQmRsDmfnFcDE3hy5h8EtEtdr78Nl/2sne1uV7o/FLWJzsgZj5y3F8vDUCA4Jc0at5QzlZw97K7K5fJyOvCHsik7HjXDJ2nU+SKx8JIovD410aY2LPJnC3Z95eIjI8Dtbmcn7DxJ5+WLEvBsv3Rsl0aC+uOY5Z607Ka/CgFm5y8rJzA4u7fh2RY1iMMz5w8RrWh1/G5fQ8+biYX3F/a0+M6+KFuOP77ulaX1cY/BJRrRO9rXMfbIkZA5vjx0Ox+G5fDBIy8+W+2ERPsRgnLC7C7Ro7oKGtBVwaWMDR2vyGMbliQQqRS1hkbBCZJS4mZ2PX+WR5wa2Y1kek/Bnf3RdPdfOFo41u9TAQESlBBJ0vDGyOp3v6ymFpPx2OkwFqyJlEuYlREB18HNGjmQu8Ha3h5WiFRo5WcLezhKnJjdfhzLwi2ekgrsUHo1Nx8OI1HItLR2GFxY6cbMzxRJfGeKKrD1ztLFFUVIS447p57hn8ElGd9kCIXLqiB2JPZAr2RiZjz4UURCXn4GhsutxuZmdpKi+iOYUlMvCtariYWHq5X0BD9AtwRUdfJ05kIyKqhK2lmczRLq7FYmhYyJlEbD+dKPfFHA2xVSQ6J0QAXKrRyKBXXIurIsqJO3libocYV6wvCwYx+CWiOicuiOI2m9gE0fsgAuHdkSm4kJiNazmFSM0pgOjMFePDxFZ+kTI2gqutBdzsLeFhb4muTZzR198VjZ05lpeI6E7GBLf0tJfbjIH+iE/LxV9ntEFwfFqevC5fSc+Ty9iXDWGoSKRTc7QxQztvRxnwimuxr7O1Xi4UxOCXiBQZGyzS84itjEjJI26rXcspwLXsQrmkspudJZxtzGHMTA1ERLXKy9Ea43v43fCYGE6WnFUgg1/R8SCGTojN1tL0hqEQ+o7BLxHpBBHgivG6YmvmqnRtiIgMjxzyYG+p+gnD6gnjiYiIiIhug8EvERERERkMBr9EREREZDAY/BIRERGRwWDwS0REREQGg8EvERERERkMBr9EREREZDCY5/c2NNfXVs3MzKyP8yHXws7NzZWvZ2ZmBjVSexvZPv3Hc6jf1H7+DKGNbJ/+K6rn92hZnFYWt1WHwe9tZGVlyX+9vb1r49wQERERUR3Gbfb29tWWMdLUJEQ2YKWlpbhy5QpsbW3rZf1q8c1FBNpxcXGws7ODGqm9jWyf/uM51G9qP3+G0Ea2T/9l1vN7VISzIvD19PSEsXH1o3rZ83sb4gB6eXmhvok3ihovaIbURrZP//Ec6je1nz9DaCPbp//s6vE9erse3zKc8EZEREREBoPBLxEREREZDAa/OsbCwgJz5syR/6qV2tvI9uk/nkP9pvbzZwhtZPv0n4UOv0c54Y2IiIiIDAZ7fomIiIjIYDD4JSIiIiKDweCXiIiIiAwGg18iIiIiMhgMfnWYr6+vXFWu4vb6669DjQoKCtC2bVvZxmPHjkEtHnzwQTRu3BiWlpbw8PDAuHHj5IqBahETE4OJEyfCz88PVlZWaNq0qZzdW1hYCLX44IMP0L17d1hbW8PBwQH6bvHixfJ8ifdkhw4dsGfPHqjJ7t278cADD8hVnsT1ZMOGDVCLDz/8EJ06dZIrjrq6umLEiBGIiIiAmixZsgStW7cuXxihW7du+PPPP6FW4pyK9+mMGTOgBnPnzr0lbnF3d4euYfCr4959911cvXq1fHvzzTehRq+++qr8Y6U2/fr1wy+//CL/QK1duxYXL17EI488ArU4d+6cXAL8q6++wunTp7Fo0SIsXboUb7zxBtRCBPKjR4/GlClToO/WrFkj/8jOnj0b4eHh6NWrF4YOHYrY2FioRU5ODtq0aYMvvvgCarNr1y4899xzOHjwIEJCQlBcXIzBgwfLNquFWFH1o48+QlhYmNz69++Phx56SF5f1CY0NBRff/21DPbVpGXLljfELSdPnoTO0ZDO8vHx0SxatEijdlu2bNEEBgZqTp8+rRFvyfDwcI1abdy4UWNkZKQpLCzUqNXHH3+s8fPz06jNd999p7G3t9fos86dO2smT558w2Pis/f6669r1EhcT9avX69Rq6SkJNnGXbt2adTM0dFRs2zZMo2aZGVlaZo3b64JCQnR9OnTR/PCCy9o1GDOnDmaNm3aaHQde3513Pz58+Hs7CyHBIjbr2q6nSwkJiZi0qRJWLlypbytrGapqan48ccf5S10MzMzqFVGRgacnJyUrgbdRFw7jhw5InsKKxI/79+/n8dLTz9rglo/byUlJfj5559lz7YY/qAmogd/+PDhGDhwINQmMjJS3skVw6see+wxREVFQdeYKl0BqtoLL7yA9u3bw9HREYcPH8asWbMQHR2NZcuWqeKwiY6Z8ePHY/LkyejYsaMcP6pGr732mrwFm5ubi65du2LTpk1QKzGs4/PPP8cnn3yidFXoJikpKTKYcHNzu+Fx8XNCQgKPlx5eP2fOnImePXsiODgYaiJuk4tgNz8/Hw0aNMD69evRokULqIUI6I8ePSqHPahNly5d8MMPP8Df3192br3//vuyw0cMWxEdebqCPb86MBj85k2McxJefPFF9OnTR44HeuaZZ+RYyuXLl+PatWtQQxtFkJSZmSmDerWeQ+GVV16R4yu3b98OExMTPPnkk/IPl5raKIiJfEOGDJHjY8X7VW3tUwvRtorEe/Hmx0j3TZs2DSdOnMBPP/0EtQkICJATn8XYZjHW/qmnnsKZM2egBnFxcbJja9WqVXLSqdoMHToUo0aNQqtWrWSv9ubNm+Xj33//PXQJlzdWoPdFbLfL8lDZh+Ly5ctyMoC4IIhvV/reRnE75I8//rjhD6/omRIB4uOPP65zH5baOIfx8fHw9vaWt5l1+TbenbZRBL5icp94X65YsQLGxrr9vfpuzqFol5gslp6eDn0d9iCGFv366694+OGHyx8Xf4hFoCEmU6mNuLaIXkORFUFNnn/+eZnFQmS2ELeW1U4EUSKTjJhYq+/EeROfP/F3ruLfPfFeFddNkfmo4v+pwaBBg9CsWTOZyUNXcNhDPXNxcZHb3RC9h4JImaWGNn722WfylkgZEUDdd999cka6Lgf393IOy3p8xQVOl91JG8WXMhH4irRZ3333nc4Hvvd6DvWVubm5PEciS0DF4Ff8LGbTk+4T1w8R+IqAfufOnQYR+Ja1W9evmTU1YMCAW7IfTJgwAYGBgXKInNoC34KCApw9e1ZmltElDH511IEDB2QPrwgq7O3t5dggMQyiLG+sGtzcDjG2SxDf8EUPt74T47TFJsbkiXHbYtD/22+/Lduny72+d0J8Yenbt688lwsWLEBycnL5/+libse7IdKAicmK4l/RQ1OWh1r0ZJS9Z/WFGCMqck2LMfbiPSjSLIl2iXH3apGdnY0LFy6U/yzmSYhzJiaF6fu1U0ySWr16NTZu3Chz/ZaN1RZ/I0SebTUQaRLFrXNxhywrK0uOjxWB/tatW6EG4rzdPEbbxsZGjodVw9jtl19+WebZFp+1pKQk2cElhjeKoSs6Rel0E1S5I0eOaLp06SJTK1laWmoCAgJkCpGcnBzVHrLo6GhVpTo7ceKEpl+/fhonJyeNhYWFxtfXV6aZio+P16gp/Zc4Z5VtavHUU09V2r4dO3Zo9NGXX34p0yiam5tr2rdvr7o0WeK8VHa+xHnUd1V91sTnUC2efvrp8vdnw4YNNQMGDNBs375do2ZqSnU2ZswYjYeHh8bMzEzj6empGTlypExjqms45peIiIiIDIbuD84jIiIiIqolDH6JiIiIyGAw+CUiIiIig8Hgl4iIiIgMBoNfIiIiIjIYDH6JiIiIyGAw+CUiIiIig8Hgl4iIiIgMBoNfIiIiIjIYDH6JiIiIyGAw+CUiMgDJyclwd3fHvHnzyh87dOgQzM3NsX37dkXrRkRUn4w0Go2mXl+RiIgUsWXLFowYMQL79+9HYGAg2rVrh+HDh+PTTz/lGSEig8Hgl4jIgDz33HP466+/0KlTJxw/fhyhoaGwtLRUulpERPWGwS8RkQHJy8tDcHAw4uLiEBYWhtatWytdJSKiesUxv0REBiQqKgpXrlxBaWkpLl26pHR1iIjqHXt+iYgMRGFhITp37oy2bdvKMb8LFy7EyZMn4ebmpnTViIjqDYNfIiID8corr+C3336TY30bNGiAfv36wdbWFps2bVK6akRE9YbDHoiIDMDOnTtlVoeVK1fCzs4OxsbGcn/v3r1YsmSJ0tUjIqo37PklIiIiIoPBnl8iIiIiMhgMfomIiIjIYDD4JSIiIiKDweCXiIiIiAwGg18iIiIiMhgMfomIiIjIYDD4JSIiIiKDweCXiIiIiAwGg18iIiIiMhgMfomIiIjIYDD4JSIiIiKDweCXiIiIiGAo/h87uEu1AdCFoAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "x = np.linspace(-5, 5, 100)\n",
    "ax = plt.figure(1, (8,3)).gca()\n",
    "plt.plot(x,  np.sin(x)-0.25, x, 0*x,'r--')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('f(x)')\n",
    "ax.set_xticks(np.arange(-5, 6, 1))\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because `plt.plot` does not by default draw a line for $y=0$, a dashed, red line is superimposed to provide more obvious intersections between the function $f(x)$ and the horizontal axis $y = 0$.  It appears (from inspections) that the roots of interest are near -3.5, 0.3, and 2.8.  One could certainly zoom in somewhat to improve the precision of the estimates, but arbitrary precision will require other techniques.\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**:  A simple, somewhat realistic model of a spherical nuclear reactor consists of two  regions: (1) an inner spherical volume of *fuel* of radius $R$, in which neutrons induce fission, leading to energy release and additional neutrons, and (2) an outer, infinitely large *reflector* surrounding the fuel, which minimizes the number of neutrons that escape the fuel.  The reactor is *critical* if the number of neutrons produced and lost is constant over time.  For this simple, two region model, the condition for criticality is\n",
    "\n",
    "$$\n",
    "  BR \\cot(BR) = 1 - \\frac{D_r}{D_f} \\left ( \\frac{R}{L_r} + 1 \\right ) \\, ,\n",
    "$$\n",
    "\n",
    "where $D_f$, $D_r$, and $L_r$ are fixed, material properties of the fuel and reflector.  Hence, in order for this criticality condition to be satisfied, an appropriate value of $B$ must be determined.  Here, $B$ is called the *buckling* and is a function of several fuel material properties.  For $R=20$, $D_f = 0.9$, $D_r = 1.1$, and $L_r =  1.7$, estimate the first positive value of $B$ that satisfies the condition.\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**: The solution of transient heat-conduction problems often involves a technique called *separation of variables*.  Application of this technique to an infinite slab of thickness $2L$ subject to certain initial and boundary conditions leads to the transcendental equation\n",
    "\n",
    "$$\n",
    "  \\cot \\lambda L = \\frac{\\lambda L}{hL/k} = \\frac{\\lambda L}{\\text{Bi}} \\, ,\n",
    "$$\n",
    "\n",
    "where $\\lambda$ is a dimensionless, undetermined parameter, $h$ is the heat transfer coefficient, $k$ is the thermal conductivity, and $\\text{Bi}$ is the *Biot number*.  The Biot number quantifies how hard it is for heat to flow *within* a body relative to how hard it is to flow *through the outer surface* of the body.  Only for certain values of $\\lambda$ can the equation be solved. Determine these values graphically by plotting the left-hand and right-hand sides of the equation as functions of $\\lambda L$.\n",
    "\n",
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determining Roots Numerically\n",
    "\n",
    "The numerical solution of $f(x) = 0$ can be challenging since it requires we know at least a little about the desired solution, particularly the range in which the root is contained. We'll explore three basic schemes for solving $f(x)$: the [bisection method](http://mathworld.wolfram.com/Bisection.html), [Newton's method](http://mathworld.wolfram.com/NewtonsMethod.html), and the secant method.\n",
    "\n",
    "### Bisection\n",
    "\n",
    "Application of bisection to finding a root of $f(x)$ requires a key piece of information: a range $x \\in [L, R]$ in which a *single* root lives.  The algorithm is actually quite similar to [binary search](A12_Searching.ipynb).   A central point $C$ is selected, and the function is evaluated at that point.  Then, if $f(C)$ has the same *sign* (i.e., $+$ or $-$) as $f(A)$, we know that the root must be between $C$ and $B$.  Why?  If there is only one root $x^*$ between $L$ and $R$, then $f(x)$ must have a constant sign for $x \\in [A, x^*)$  and the opposite sign for $x \\in (x^*, B]$.  If, instead, $f(x)$ changes signs between $x=A$ and $x=x^*$ then it *must* cross the $x$ axis, at which point there *must* be a second root.  \n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**: Write a `lambda` function for $x$ that returns `True` if $x \\geq 0$ and `False` otherwise.\n",
    "\n",
    "<hr />\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An important difference between bisection and binary search is that now the search can require an unlimited number of comparisons because the search range $[A, B]$ is split in half repeatedly until a root is *converged* to within a desired tolerance $\\tau$.  In other words, roots will very rarely be found *exactly*; rather, we often must accept *approximate* values, as was first introduced [earlier](A04_Loops_in_Python.ipynb). \n",
    "\n",
    "In pseudocode, the bisection method can be written as follows (where the notation from binary search is adopted as closely as is reasonable):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "'''Bisection method for finding an isolated root of f between a and b'''\n",
    "Input: f, a, b, tau\n",
    "\n",
    "# Define the left and right boundaries\n",
    "Set L = a\n",
    "Set R = b\n",
    "\n",
    "# Define the central point\n",
    "Set C = (a+b)/2\n",
    "\n",
    "# Go until f(C) is close enough to zero (i.e., that\n",
    "# C is close enough to the root between a and b\n",
    "While |f(C)| > tau \n",
    "    If sign(f(C)) == sign(f(L)) then\n",
    "        # the root must be between C and R\n",
    "        L = C\n",
    "    Otherwise\n",
    "        # the root must be between L and C\n",
    "        R = C\n",
    "    C = (L+R)/2\n",
    "Output: C\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Application of this algorithm to find the first positive root of $\\sin(x) - 2/5$ is shown graphically below for the first two iterations.  At each step, the search space is reduced by half, illustrated by the shrinking highlighted rectangles.  By the second iteration, $C = 0.375$, which is visually quite close to the root of interest $x = \\sin^{-1}(2/5) \\approx 0.4115$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T18:55:46.721577Z",
     "iopub.status.busy": "2026-07-05T18:55:46.721386Z",
     "iopub.status.idle": "2026-07-05T18:55:46.843984Z",
     "shell.execute_reply": "2026-07-05T18:55:46.843440Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAG5CAYAAACTEQDKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZIVJREFUeJzt3QV0k1cfBvCn3uLuxYu7u/sYYwwbUNwpUNyluDst7g7DBgzZcPfh7u4O1Xznf1n6taWStmljz++cQJPevLl9kzRPr1ppNBoNiIiIiChc1uF/m4iIiIgYmoiIiIh0xJYmIiIiIh0wNBERERHpgKGJiIiISAcMTUREREQ6YGgiIiIi0oGtLoUoYgEBAXjy5Anix48PKysrnjIiIiITIMtVfvz4EWnSpIG1dfhtSQxNeiKBydnZWV+HIyIiolj08OFDpEuXLtwyDE16Ii1M2pOeIEECfR2WiIjIbDx8CLx+HfuPmzQpEFa7xocPH1Sjh/ZzPDwMTXqi7ZKTwMTQREREFNyDB0CRIsC3b4h1jo7A9etA+vRhl9FlaA0HghMREVGMe/XKMIFJyOPK40cXQxMRERGRDtg9Z4TNl/pIw5GVLFn4zZZERESWjqHJyAJT9uzG299LRERkyRiaDMTX1xf+/v7Bbnv82ArfvjkYpD4S1B4/9kaKFBqDPD6RsbKzs4ONjY2hq0FERoChKZbJ1MZXr17B29v7h+89eeIIIFNsVynI4z9BokQGGqVHZKRkRk3ChAmRKlUqLlxLZOEYmmI5MD1+/Bjx4sVDsmTJ1F+wQac4vntn2JXEZTXUTJnY0kQUdKXgz58/4+XLl3ByckKiRIl4cogsGENTLJIWJglMsuJoaOtBOBimZy7I4zuosU1E9H8SlqRl+MWLF6rFidskEVkuLjkQi2OY5Bcvf+kSmR5ZsFbGIIYch0hEloWhKZZof9lKlxwRmRZb2++N8n5+foauChEZEENTLGPTPpHp4fuWiARDExEREZEOGJqIiIiIdMDZcxRIo/l+iYolS5agdetWuHPnLjJmzKjz/e7du4fMmTNh0aLFaNmypbrt6NGj2L17N9zd3X+Y4l2xYgX1/759+yNdx+jclyyb9r0hQxM5FpwoaszhvcPQRMFe0FEd51q9+k84dOgYkidPHaljSHm5X+bMWQLvd/jwUYwY4YFmzVoiXrzgoWn6dE/1f1TqqQ2EHMtLkSWvmYAAWWsN+PqV548oKuT9Y+oYmiiYUJaP0kmKFMnVJbIcHR1QokSJUOsg/4esT+7cuaJWwVCOTxSZ14xcZDeV/ybSEVEkmcN7h2Oa6IcPhqhcli1bAjs7K9y/f09dr1y5AgoUyIPTp0+hQoWySJAgDrJly4wJE8ZBowkIvJ+Ul/vJ/eX6iBHD0a9fH1UfF5dM6ntyOXBgf+Bx5RL0sUeO9ECpUsWRIkUSJEmSAMWKFcLixQulbSlYOX38nLxY7jkQEpp44TngawBROgfWZpA4zCD3kbF69uwZmjdvih49emHIkGHYsmUTBg0aoLZrcXVtHup92rRpi7dv32DWrJnYsGEjUqVKrW7PlSvsFiYJXu3bd4Czc3p1/cSJ4+jevavasmbIkKEx9NMREZGlYWiiGPP69Wv8+ecOFCtWTF2vUqWKajFavXpVmKFJtpjRhp8CBQrqNKh84cLFgV8HBASgQoUKas+wmTOnY/DgIVxjh4iI9IKhiWKM7AqvDUxaefPmw7//ntfr4+zduxfjxo1RXYGyKXJQsl9YypQp9fp4RERkmcygh5GMVdKkSUPdFPirHqcfnTx5EjVrVlNfz507HwcPHsHx46cwYMAgdZs+H4uIiKK+/+r27Zth6tjSRCZt7do1aj+/rVu3wdHRMfD2LVtM/81JRGTq/Pz8sGbNckyYMAr37iUEUBemjKGJjI60RunaSiR7gslmqjYyNeM/cr+VK5fHaB2JiCj8TerXrVuFceM8cPfubXVbokSV8e4dTBq758jo5MmTV/0vA7mPHTuG06dP4+PHj6GWrVXrJ3z69AnNmjXBnj17VMuTLHGgDV5ERBR7AgICsHHjOhQvngcdOjRXgSlZsuQYNWoi1q790+SfCoYmMjoy+61fvwHYtu1PlC9fBiVKFMWZM2dCLVupUiUsWLAIly5dRN26P2PIkEGoV68++vbtH+v1JiKyVBqNBtu3b0GZMgXRsmUj3LhxDYkSJcbw4WNx4cIddOvWG05OTjB1Vhr5SSnaZNZWwoQJ8f79eyRIkOCH73/79g13795FpkyZgo29CersWaBwYcM9GSdPAoUKGe7xiYyVvH/v3buLVKkywd4+9PcvkSXSaDTYt+9vjBw5GGfOnFS3yWegm1svdO7sHuzz8Px5oFw5w9VV/vYO7TMuos/voDimiYiIiCLt+PEjGDFiEA4fPqCux4kTB506dUfXrr2RJEkSszyjZts95+npGdiqU7hwYRw6dEin+x05ckQNLC5QoECM15GIiMjUXLz4Lxo0qI1q1cqowGRvb6/CknTDDRs2xmwDk9mGprVr18Ld3R2DBg3CuXPnULZsWdSsWRMPHjwI937SNNe8eXNUrlw51upKRERkCm7duolWrX5H6dIFsGvXdjVruUWLtjh//hbGj5+GFCnMfyFhswxNU6ZMQZs2bdC2bVvkzJkT06ZNg7OzM7y8vMK9X4cOHdCkSROULFkywsfw9vZW/aBBL0RERObm6dMncHfviKJFc+KPP9ao2377rTFOnryCmTPnI106Z52OI+sdG2piswwlTpYs+scxuzFNPj4+aqZV//7BZ09Vq1YNR48eDfN+ixcvxu3bt7FixQqMGjUqwscZO3YsPDw89FJnIiIiY/P27VtMnz4BXl7TA9fNq1atFoYOHY18+SI/hMXZ+fuEp9evI1+XgABZKFMGmQNBluXTmQSm9N+3NY0WswtNr169UotqhdxvTK4/e/Ys1PvcvHlThSwZ9yTjmXQxYMAA9OzZM/C6tDRJa1bs0qAgzqEoTiE3LsMJX/EVTriM3DiFojiHgjJBMpbrREREpkwC0rx5szB58li8e/dW3Va8eCl4eIxDqVJlo3Vs+ZiMykelv//30JQ4cdRCk76YXWgKulJ0yGmRIW8TErCkS05ajbJly6bz8WXxREMtoGgLX7TDfLhjGrLhJgJgBT/YwgoaaGAFW/jBGhrcgAumwR3z0Q5+sDNIXYmIyDTI5+Hq1cswevRQPH78SN2WM2duNbi7Zs2fQ/0MtTRmF5qSJUumBqeFbFUKa7d7WWlaVpyWAeNubm6BK5pKyJJWp927d6sFFGOn7t/7Xb99C7tMblzCCjRDPlwIvE0Ckj18fyibFbcwC25oj3lohhW4jDwx3t9LRESmRT7vdu7cjuHD++Pq1cvqNhmnNGjQCDRu7BpsmypLZ3ahSaY+yhIDsqXGr7/+Gni7XP/ll19+KC8LWV28ePGH5Qr27t2LDRs2qGULYov0t16/Ll2MwW+X5UelaTL+2QPI7l4TVn4+sPKPeE1SCVMiv80l/GtbDDen/4WPhcrHaH8vERGZjtOnT2LIkD44cuSgui6rePfuPRDt27uFuRCzJTO70CRkrJGrqyuKFCmiZsLNmzdPLTfQsWPHwPFIjx8/xrJly2BtbY08eYK3wKRIkUK9WELeHhskuIQMLxKa/M5fgm2PmoCvN6xkRFwkWPn7w1rjjew9aiLg2EnZ3E2/lSYiIpNy585ttTDlxo1r1XUZbiJrLfXo0R+JZeAQWc6SA40aNVLLDIwYMUItUnnw4EHs2LEDGTJkUN9/+vRphGs2GRVfX9i2bAb4+kQ6MGmp+/n6wrqFHOfHrrzYJF2hFSuWR5IkCWFra4Xp06ep58fDY3is12XWrJnInj0rnJzsVV3evXuH1q1bIkuWjHp7jCtXrqif7d69ezqVX7p0iapL0PKrV69S58nQwquH1NkQzyER6e7169fo189dLR8ggUnGKTVp0gJnz97AiBHjGZgiwL3njGjvubBoZnsCXd3kyYp2PTVWVtDMmAVNp84wlCJFCuLz58+YMmW6eoNmzJgRY8aMgqfnbPj5xd5WiOfPn1d1adOmLVxdW6gxbEWLFlVhRZ7PggVl9mH0/fHHBjRq1AB//71PbUYckZcvX6rlL+TxtZMN6tSpjcuXL+H2bd2CV0wJrx7Hjx9HunTp1MXccO85MofXsMyImzhxlPqcEpUrV1dBKW/e/DB2/jE4e457z5kTCUp6bmGwmjkdmo6dZIohDOHSpUto27adWqU9pn358kXthxSaK1e+D3hs06YdihUrFnh7lixZdBo4Kb+EYmLX7uTJk6tLbE0t1tfPUKJECb0ch4j0R35X/fHHWnh4DMD9+9//2JGQNHLkRFSqVJWnOpLMsnvOrJw7B6ubN/XSyiTkOFY3bnzfblqPbt26hTZtWiFHDhfEjx8H6dOnxS+//BxskL2228nPzw9z5nipr+Ui3WHSyiS0twXtnpI3vZeXJwoXLoB48ZyQLFliNGxYH3fu3AlWh0qVKiB//jyqO7ZMmVKqHm3btg61vlK2efNm6utSpYoH1kOE1j0n3+/WzQ1z585Bnjw5ESeOA5YtW6q+Jz9LoUL5kTBhPCRKFB+5c+fAoEEDA39maWUSVapUDPzZ5HZdu+ekrjt2bMf9+/eDnZ+gC7qOHj1KPa7UK1Wq5Oq5kBaroORnkpaiTZs2qha2uHEdMWLE9wVa5fxXqFAOqVOnQIIEcVGgQF5MnDgBvkG6ciOqR2jdcxKQf/31F/WcyePJc6g9b1r79+9X912zZjUGDx4EZ+c0SJw4AapVq4LrMjOCiKLkxImjqFy5JFq3/l0FptSp08DLazEOHjzDwBRFZjkQ3KycOqW61PQVmoQ63ulT0Oip+0k8efIESZIkxZgx45AsWXK8fftGfThKIDl9+hyyZ8+OWrV+wuHDx1CmTEn89lt99OjRS91XWlWku066seT7WqlTp1b/d+rUQQUJN7duGDNmvDr2qFEjULZsKZw9+2+wpSRkvFqLFs3Qu3dfjBo1Rg30D82sWZ7qQ1q6BRcuXIzs2XNE2LqzZctmHD58CIMHD0XKlKnUhIG1a9fAza0z3Ny6Yvz4Serxbt++hatXr6j7yM8s9Rg8eCBmzpyNggUL6dyaFbSuHTu2VwM3N2zYFOx7sjyGhBKpl/zMJUuWwoMH9+HhMQyVK1fAiROng7UknTt3FteuXcXAgYORMWMmxI0bV90u3YG//95E3SYzUC9c+Bdjx47G9evXsGDBogjrERoJPPIcyXmaNm2Gen2sWrVChdLnz5+jT5++wcrLOSpVqjTmzl2Ajx8/YMCAfqhb92dcunSVU56JIuHu3Ttq+YBNm9ar6/I+lwHebm49w2x5J90wNBm7y5flT3j9Dt6W4/3XNaUv5cqVU5egi6RJYMiXLzfmzZuLyZOnBOt2kqATtDtHG3xCdvHIOJkFC+Zj4sTJ6NHj/yuwlylTFjlzZsPUqVMwbtz4wNvfvHmDNWvWR7i2Vq5cuQKDS+7cedRMy4h8+vQJ589fDDZQct68OUiUKJEKBVpBN3yWn9fFxUV9nTNnrih1YUld5TFkfFPI+69fvw67du3E+vV/4Ndf6wXeni9ffpQoUVSFzY7SFRtkvbKLF6/8sJCrPD9Bg5hscp00aVLVYiXnXn7m8OoRmhEjhqtWMBnLpV0tv1atWmqw/ciRHmjfvoMaB6gl52fZshWB12VtmMaNG+LUqVPs+iPSgYxVmjRptNr2RN578kdc8+Zt1HpL8oceRR+754yd7Pejx1YmRY735YteDyldbmPHjkHevLnUTDQHB1v1v2xRIy0bUbV9+zY1u6Np02bqMbSXVKlSqWBw4MD+YOXlwz2mFiOtWLHSDzNLihYtpkJA06a/Y+vWLWobn9gk50eCTO3aPwc7PzJrVM5RyPOTL1++UFe+lxmNdevWQYoUSWFvbwNHRzu0bNlchd8b0p0bBfv27UWlSpV/2F6oRYuWaqzZsWP/b1UUP/9cJ9j1vHnzqf+l5YyIwibv+YUL56BAgayYPn2iCkwVK1bB4cPnMGPGPAYmPWJLk7GTrhV9D9iW4+m5ibZ3755qXEyfPv1Qrlx5FS7kr5wOHdoGbvQYFS9ePFdjmtKk+XE1d5E5c+Zg17VdejEhtGM3a+aqfmFJa1iDBr+pVpoiRYpixIhRqFo15gdZSjeXhDYJqKEJGeJSpfrxZ5DlNypUKKu6UKdOnY4MGTKqGZ6nTp1E165dovz8ydTm0M6ZjKsQb94E37VTWraC0s4cjM7rh8jc7d27BwMH9sSVK5fUdReX7Bg9ejKqV6/FbU9iAEOTscud+/s8S32S4+XKrddDrly5Aq6uzTF69JgfPrQTJkwU5eMmTZpMvfH37z8U6l5/IW+Lyb2Rwjp2y5at1EXGZckgdBlP9MsvtXH16o3AtcFictsgCRvbt+8M9fvx48eP8GeQsVpS9/XrNwar77//Rm+ygNRLxpiF9PTpk/++z317iKLq5s0bGDy4N/7668/AlbwHDvRAmzYdYWfHvUZjCkOTsStaVK+DwIUcL6BIUf0e08rqhwCzfft2tfJ6lixZI7y/vf3/WxWCDlz+6afamDBhHJ48eYwGDRrCmMlgS1lGwdfXB/Xq1cXly5dVCAn6s0WVnNvQ7i/nRwajSzda8eLFo3RsbZAK+vxJ6560nulaj9BI19zmzZvUJIE0ab63Lonly5epwahcooAo8qRlWdZamjNnhprdKuvLtWvXBf36DUWSJEl4SmMYQ5OxK1gQGhlIfOuW3ha3hByvQAHok3x4y6BjmYUmY1HOnj2DyZMn6rzQYd68edX/EyaMR40aNdUgYBl/U7p0abRr114NSJaNlcuWLafCybNnT3HkyGHkyZM32EDn2NahQzsV8mTWl3R9PX/+DOPHj1UDnGWhTKHdjmfBgnmq5Ue6vmSR05DdUeGRn1OWCvi+vEFh1fUpg9cbNWqMVatW4uefa6Fr1+5qjJX8lSk7lO/fvw916vyCunX/vwdjaKpUqapmzMm4LJmB5+39TT3Ou3dvda5HaIYMGabGXMlSCzLjMHHiJFi9eqVatmDcuAnBBoETUfjkD6NlyxZi5MjBePXq+3Ii1arVwpgxk5EtWw6evljC0GTsJOR0d1crguuLpmt3vY+TkrEw8mEtgUFmmcnUeunuGTp0sE73l+nuR48ewZw5nmo5AWnpuHXrrlot3MtrLooXL6Fm4cn3ZdyQtFxIUJGQYEgyi0/Cosxie/v2reouK126DBYvXhY4U1AC0pQp0zBz5nS1DID88pNlDmRAtK66deuuFuOUafkyQ0bOj6yeLuFy8+atmDFjOlauXK7Ov/zlKWFVxpZJyIlIjhw5sG7dHxg2bDAaNKinwlzjxk3g7t4TtWvX1KkeoZExUocOHVVltWOjcubMGemfncjSHT16CH37dsOFC9+7zCUkjRkzBdWqxfwCwRQct1ExhW1UfHxVNx0uX1Kb70aVRpYayJ0bASdOf192gIh0wm1UyBAePXqIIUP64o8/1qjr8hnTv/9wtG/fxeLGLfkbyTYqXHLAFNjZwW/JChn4A00YizVGRN3Pzg4BS1cwMBERGXlInzBhFIoUyaECk4w7bN26A86du4kuXdwtLjAZEzY3mIo8eeC/7S/Y1K4Jja8vrCIxo061MElg2vaXOg4RERkf6e7esWMrBgzogXv37qrbSpYsgwkTZiB/fv3t4EBRx5YmU1KuPAKOnVRdbDKgWw3qDkdgGemSk/uVLx9rVSUiIt3duHEN9erVxO+/11WBKU2atFi4cBV27jzIwGREGJpMTZ48CDh+CpoZs4CsWf8fjuzsVIuS+l8bplxcVDkpzxYmIiLj8/HjRzVuqWTJfPjnn11qJmuvXgNw+vQ1NGjwOxeoNDLsnjNFEow6dYZGptqfP68231V7ycnWKLLSd67c39dhkmUFYnCxRyIiinpX3IYNa9QCldoFX6tX/wnjxk3TaW07MgyGJlMmgUjWcSrIvm4iIlMhW5707u2Gw4cPqOuZMmVRYalmzdqGrhpFgKGJiIgoFsjU9rFjh6vVvGW9NlkYt1evgejWrXekl6Ihw2BoMiHeAQFwsLZiHzcRkYl1xa1fvxqDBvVSuwaIOnXqqQUq06eP2f0pSb8YmkzEQ+9vKPnvGaR3cIRHxkyoligxwxMRkZG7du0KevbsHNgVlzlzVkycOBNVq9YwdNUoCjh7zkS88vXFC19fnP30ET9duoCS589i19s36i8YIiIyLrKd1NCh/VCqVH4VmKQrbsiQUThx4hIDkwljS5OJCfjv/3P/haci8eKz5YmIyEjIH7Lbtm1Gv37d1TYoolatOhg/fjoyZMho6OpRNLGlyUT5hwhPptTydO7cOVSsWB5JkiSEra0Vpk+fhh07dsDDY7ihq2Y2smTJiNatI94UV85/0PN+5coVdf3evXswpPDqIT+X/HxExubu3Tto0KA2mjatpwKThKS1a7dizZotDExmgqHJxJlieGrXrjWePXuKlSvX4PDhY2jUqDF27tyBkSM9DF01iyPnv02btoHXr169op4HQ4em8OoxaNAQbNiwySD1IgqNt7c3Jk4cjeLFc2P37h1qb7jevQfhxInLqFnzZ540M8LuOTMNT8bcbXfp0iW0bdsONWvWjPHH+vLlC+LIgp8UqhIlSsTKmdHn85AlSxa9HIdIHw4e3IcePTrh5s3r6nr58pUwZYonXFyy8wSbIbY0mRlDtTzdunULbdq0Qo4cLogfPw7Sp0+LX375GRcvXgwss3TpEtUd5OfnhzlzvNTXcpHuFk/P2aqM9ja5aFsZpO5eXp4oXLgA4sVzQrJkidGwYX3cuXMnWB0qVaqA/Pnz4ODBgyhTppSqR9u2rcOss9y/SZPGcHZOgzhxHJAmTUpUrVoZ58+fD1Zu3bq1KF26JBIkiIuECeOhZs3qqosxpBMnTqifOUWKpIgb1xHZsmVBz57uwcocPnxYPUaiRPFV/aSe27dvD1ZGe5727duHLl06IWXKZOqY9evXw5Mn31cO1vL19UW/fn2RNm0qdbxy5crg5MmT0FXQ7jl53EaNGqivq1SpGPg8yO1af//9t6p/4sQJ1OOVLVsa//zzT7BjyvHkfmfPnlXPkzxfci7E6dOn1TmX7jV5LuX/pk1/x/3794P9/OHVI7TuOdkVfuDAAciaNROcnOzV669r1y549+5dsHJyvzp1amPnzp0oWrSQqkPu3DmwePEinc8ZkXj58gXat2+O2rUrqcCUPHkKLFiwElu3/s3AZMYYmiwoPP399k2MPZ58mCdJkhRjxozD9u07MXPmbNja2qJUqeK4fv37X2C1av2kuoPEb7/VV1/LRbpb5LrQ3iaX1KlTq9s6deqgwkelSlXwxx+bMWuWJ65cuYyyZUvh+fPnwerx9OlTtGjRDL//3gR//rkDHTt2DrPOP/9cC2fPnsG4cROwc+cezJ7thQIFCgb7oB07doz6UM+VKxdWr16HpUuXq72iKlQoq8bdaO3atUvd9vDhA0yaNAXbtv2FgQMHB6vfgQMHULVqJXz48B7z5y/EypWrET9+fNSt+7MKZiF16NBWNfMvX75K1fHAgf1o3rxZiDLtMGXKJDRr1hwbN25BvXq/oUGDenj79m2kn0N5fkaNGqO+ludP+zzI7WLlyhWoWbMaEiRIgMWLl6rzkSRJEtSqVf2H4CSkHrIdxJo16+HpOUfddv/+PWTLlh1Tpsg4tl0YO3a8es5KlCiKV69e6VSPkCRU16tXV52Hpk1dsXXrdri798SyZUvV+Zauk6AuXPgXffv2QvfuPdQ5y5s3H9q1a6PCNlFEAgICsHTpAhQunANr1ixXLflt23bGmTPX0bBhE6Nr2Sf9YvechYSn058+osftW7hYpFiMPE65cuXUJfBx/f3Vh1y+fLkxb95cTJ48BcmTJ1cXkTJlymBdQ3I9tO6i48ePY8GC+Zg4cTJ69OgZeHuZMmWRM2c2TJ06BePGjQ+8/c2bN+pDulKlSuHW9/Xr1yrMyYd306b/DyK//lov8OuHDx/Cw2MYunRxw7RpMwJvr1KlqmpRkzE3q1d/DzvdunVB+vTpcfToiWAr+7Zs2Srw64ED+yNx4sT455/9iBcvnrrtp59qqxa0vn17o0GDhsF+4VavXiPY48rP1r9/Xzx79gypUqXCtWvXVDCQD//x4yeoMlWrVlXn0tW1KSJLnhsXFxf1dc6cuYI9F9K91qNHd1XfP/74/3iiWrVqqRabwYMHonLlE8GO5+raAsOHBx+nJuFYG5C1rxM5prTyrV69Cl27dgu3HqHZvXs3du/epYJl7959As+Ds7Mzfv+9EZYvX6a6g7UknB08eEQ9X0Jet3v3/oM1a1YFew0ThXT16mV0794Bx48fUdfz5SuAadPmokgM/V4l48OWJjNn89//MsZpagxuAildbtIqkzdvLtU94uBgq/6/efMmrl27GuXjbt++TQUJCTbyGNqLhIZ8+fKr1pegJJREFJiEtJDI2JjJkyeq4CXdbfIXZFDyQSyPJa04QR9bQlG5cuUDH/vGjRu4ffs2WrVqE+ZWCJ8/f8bJkydUYNAGJmFjY6NaRx49ehTYIqdVu3adYNelRURou7L279+n/m/SJHhAkvAlrXz6dPToURXamjdvEexcyDmrVq0GTp8+pX7GoKTVK7S1a/r374fs2bOq14hcpMtT7hvV18m+fXvV/y1aBJ8tWL9+A8SNG1cFoqDy5y8QGJiEPGcuLtmCdRESBfX161d4eAxE6dIFVGCS19WYMZOxf/8pBiYLw5YmMw5L0spUMJYGhPfu3VONS+rTp58KFBJerK2tVReT/MKJqhcvnqvuF2mJCE3mzJmDXdd26UVEzsXu3f9g1KgRmDRpAvr06aWClASQkSNHq24zbdeadB2FRn4+8fLlS/V/unTpwnw86S6TnyNVqh/rlyZNmsDWr6CSJk0a7LqDg4P6X3s+teUlQAYlgSnkfaNLngchY5TCIqFKPkzCey6aNWuiQox0yRYpUlR19clzIV2lUX2dvHnzWv3M2lZMLTmunJuIzqv23EbndUrma+/ePXB374h7976Pofzpp18wYcIMODv/P3iT5WBoMjOxHZa0ZLyLq2tzjB79fSxK0K6QhAkTRfm4SZMmU/Xfv/9QYGgIKuRtkflZM2TIoMYWaVuL1q9fhxEjhsPHx0eNwUmWLJn63rp1G8LdH0r7YS2tRWHRhkhZaiEk7eBu7ePpSvvhL911adOmDbxdWoBCBoXokudBTJ8+E8WLh95dpu1iDeu5eP/+vWo5HDJkGPr16x94u4w5ksAVVTKWTn5mCa9Bg5OEVDk3Es6IojLQe8CAnli3bqW6niZNWkyaNAu1a9flybRg7J4zs244CUvb8+TDsQKFUD1xklgblCiPEzLAyKywx48f63R/e/vgrShaMt5FPvyePHmMIkWK/HDJmzevXuqfLVs2DBo0WB1PZn2JatWqqxYM6XoL7bHlor2vdPUtWbLoh0HHWtICU6xYcWzatDHYzyjdW6tWrVCtVHKcyChfvoL6f9Wq77/UtST8SYiIirCeh9KlSyNRokRq/aSwzoW9vX2ErxF5LkO+ThYuXKDGNulSj9BUqlQ5MLgHtXHjH6rbT/t9Il3Ia3TlyiUoUiSnCkzyuu3YsRtOnrzCwERsaTJ1hmpZCknCjUwJz549hxp7I7PSZLxQeF1WQWnDz4QJ41GjRk011idfvnzqw7pdu/ZqOQOZrl62bDkVQKTF5siRw8iTJy86duwU6fpeuHAB3bq5qXEvWbO6qA98GRsjt/ft+70VJGPGjBg+fASGDBmkVvqVgdmJEiVWXVUyrV/qoR3oPGPGbDULrnTpEmpgtjTdy0w6GRe1fPn3UDN69FjUqFEVlStXRK9evdVjylIKsm6VzKSL7POWM2dONdZrxoxpapZd5cpVcPnyJTWLTLq9oiJPnjzq/wUL5qkuShnvkylTJtWqJa1MrVq1UK1CMjZLpli/evUS//77r/pfZh+GR+okz5+8LqRVTVZLPnjwABYvXqgCma71CEkGfUvAHTCgHz58+IBSpUrj4sULahB/wYIF0ayZa5TOBVmeW7duwt29g1p7SeTNmx8zZsxH4cJsraTv2D1noowlLGlNnTpdfXCPHz9WDfYtWLAQ1q/fiKFDB+t0f1ki4OjRI5gzx1ONM5K/9m7duquCi5fXXNUlJLPw5PvSOiPjgOTDsWjRqM1akbEu0jokx5NZcnLuZHyUzNJzc+saWK5//wFquYEZM6ZjzZrVqiVJ7itdPu3bdwwsV716dezbd1DV3d29m1o3SAJj0MHc5cuXx549e9WHuaw1JD+HDGbftGkrateuHaWfQ7oXU6RIiWXLlmDWrBlqkPO6dX+gadPGUTqeBBOZUThz5nRUrlxBtQAtXLhYDbKWgCZhUMaAyTIQsvRCihQp1GOGHIQdlhUrVqlZeDILUFrD5DmU5R7q1PlJ53qEJM/dxo2b1fpQS5cuxtixo1Uok7AkSxeE1q1LFHK9sxkzJmHcOA/1HpfNdQcO9EDnzu7q9xqRlpXGmPfbMCHyF27ChAnVuI3Q/sqXD9G7d++qD4OwZliFRZ6hU+8+ovi/Z1R/qszxMuYVv4nMjbx/7927i1SpMsHePnLvXzJuZ86cgptbG1y+/H0h3ooVq2LqVC9kzsyV542Jv7+M15TxoTLrOHY/v4NiS5OJSG5nh5R2dnB2cGRYIiKKJhnvNmrUEHh5TVetvjKhYOzYqWjcuBn/EKUwMTSZiHQOjrhdtCQcrK34hiYiioZ//tmtxi7JCvWiYcOmGDduKpIlC75sBVFIDE0mxMHaGuyJIyKKGpnEMHBgT6xatVRdlzF6U6fOQbVqMb95OJkHLjlARERmTYbubtq0HkWL5lSBSbuMwIkTlxmYKFLY0hTLOO6eyPTwfWu6nj59gl69umDbts3qeo4cuTBz5gIUL17S0FUjE8SWplgi01blr5uQ+3MRkfH7+vWL+t/WltPPTSnoLl++CMWK5VKBSX4H9+8/DIcOnWVgoihjS1MskcUaZUqjbPUg64DItEZZbVqX5QJkyQGZailFOaaJKHY/eCUwyZYaceMmgrW1nuc6U4yQAd7durXHvn171PVChYpi9uyFyJ1bPzsIkOViaIpFsiiiLJr24sULtS6EriQ0BQQwMBEZigSmJEmCb4xMxkeWDpg3bzY8PAaoVn1ZE2/w4JFqkUr5I5UouvgqikXSqiTbRUiLk6xwrOv+YLKol2QsWdBL34t6EVH4pEuOLUymsQVKly6tcezYYXW9dOlyauySbJNEpC8MTQYKT/JXj65/+Uhokn1LpThDExFR0N+P/vD0nIaRIwerldvjxYsHD4/xaNOmI6ytOWyX9IuhiYiITNL161fRqZNs5n0icAuUmTPnI336DIauGpkpxnAiIjIpMrRh6tTxKFOmoApMMuRBBnpv3ryLgYliFFuaiIjIZFy9elm1Lp09e0pdr1atFqZPn4u0adMZumpkAdjSREREJtG6NHnyWJQtW0gFJmld8vJagvXrtzEwUaxhSxMRERm1a9euoGPHloGtS9Wr/6Ral9KkSWvoqpGFYUsTEREZ/dglCUyyZMucOUuxbt2fDExkEGxpIiIio3PjxjXVuqSdGSdjl2bMmMewRAbFliYiIjKqdZdmzZoSODNOtpzy9Fykxi6xO44MjS1NRERkFG7fvoXOnVsFrupdqVI1zJq1AOnSORu6akQKW5qIiMgo9owrXTq/CkyyqrcM9N60aScDExkVtjQREZHBPHz4AG5ubbBv39/qetmyFeDpuRgZMmTks0JGhy1NREQU6zQaDVauXIKSJfOqwOTk5IQJE2bgzz//YWAio8WWJiIiilUvXjxHt27tsWPHVnW9aNESaikBF5dsfCbIqLGliYiIYs3WrRtRvHgeFZjs7OwwfPhY7Np1iIGJTAJbmoiIKMa9e/cOfft2w5o1y9X1PHnyYd685ep/IlPBliYiIopRMmZJxi5JYLK2tkavXgOwb99JBiYyOWxpIiKiGPH161cMG9Yfc+bMUNczZcqCefOWoXjxUjzjZJIYmoiISO/OnTuDdu2aqe1QRJs2nTBq1ETEjRuXZ5tMFkMTERHpdZPdKVPGYdw4D/V1qlSpMXv2IlStWoNnmUweQxMREeltG5T27V1x6tRxdf3XXxtgyhQvJE2alGeYzAIHghMRUbQXqlyyZD7KlCmgAlPChAkxf/4KLFmyloGJzApbmoiIKMpevnwBN7e2+OuvPwO3QZGFKp2d0/OsktlhSxMREUWJBKUSJfKq/+3t7TF69CS1DQoDE5krtjQREVGkfP78GYMG9cKiRXPV9dy586ruOC5USebObFuaPD09kSlTJjg6OqJw4cI4dOhQmGU3btyIqlWrInny5EiQIAFKliyJXbt2xWp9iYhMwdmzp1G2bKHAwOTm1pMLVZLFMMvQtHbtWri7u2PQoEE4d+4cypYti5o1a+LBgwehlj948KAKTTt27MCZM2dQsWJF/Pzzz+q+REQE+Pv7Y+LE0ahSpSRu3bqBNGnSYuvWvzFmzGT1xymRJbDSyLQHM1O8eHEUKlQIXl5egbflzJkTdevWxdixY3U6Ru7cudGoUSMMHTo01O97e3uri9aHDx/g7OyM9+/fq9YqffL3B96+BWxtARsbvR6aiChCDx7cV0sJHD16KHApgalT5yBJkiQ8exQr5HPQzw9InFj/n4Py+S0zPnX5/Da7liYfHx/VWlStWrVgt8v1o0eP6nSMgIAAfPz4MdxfCBK+5CRrLxKYiIjMzbp1q1CqVD4VmOLFiwcvryVqKQEGJrJEZheaXr16pZqRU6ZMGex2uf7s2TOdjjF58mQ10LFhw4ZhlhkwYIBKpdrLw4cPo113IiJjIb/X2rZthrZtm6q/xIsVK4kjR/5F06YtYGVlZejqERmE2c6eC/mmll5IXd7oq1evxvDhw7FlyxakSJEizHIODg7qQkRkbo4fP6L2jbt//x5sbGzQt+8Q9OkzCLYyRoDIgpndOyBZsmTqTR6yVenFixc/tD6FNoC8TZs2WL9+PapUqRLDNSUiMi6yV9yECaMwYcJINUwhY8ZMmD9/JYoXL2noqhEZBbPrnpMF1mSJgT179gS7Xa6XKlUq3Bamli1bYtWqVfjpp59ioaZERMbj3r27qFGjnNpoVwJT48auOHz4PAMTkTm3NImePXvC1dUVRYoUUWsuzZs3Ty030LFjx8DxSI8fP8ayZcsCA1Pz5s0xffp0lChRIrCVysnJSQ3yJiIy98HePXt2UmOXZPaQzIxr0OB3Q1eLyOiYZWiSpQJev36NESNG4OnTp8iTJ49agylDhgzq+3Jb0DWb5s6dq5qlu3Tpoi5aLVq0wJIlSwzyMxARxTSZJdyrVxesWbNcXS9RorRa2TtDhow8+USWsk6TIURmnYfI4jpNRKRvp0+fRJs2TXD37m1YW1ujX7+hHOxNRsvfSNZpMsuWJiIiCp2MV5o2bQJGjRqiWthlc92FC1epViYiCh9DExGRhXj69Ila2fvAgb3qer16DTFt2lwkSpTI0FUjMgkMTUREFuCvv7ahU6eWePPmNeLEiYNJk2ahadOWXKiSKBIYmoiIzNi3b98wdGg/zJkzQ13Pn78gFi1aDReX7IauGpHJMbt1moiI6LsbN66hcuUSgYGpS5ce+PvvYwxMRFHEliYiIjMjk6JXrlyC3r3d8OXLFyRLllxttFu9ei1DV43IpDE0ERGZEZk+7e7eERs2rFbXK1SojHnzliNVqtSGrhqRyWNoIiIyE2fOnEKrVo1x794dtQfn4MEj4e7eV31NRNHH0EREZAZrL82ePRXDhvVXay+lT58BCxeu5r5xRHrG0EREZMJevXqJjh1bYvfuHer6L7/8hpkzF3DtJaIYwNBERGSiDh3aj7Ztm6pFKx0dHTFu3DS0atWeay8RxRCGJiIiE+Pv748JE0Zh/PgRqmsue/acWLJkLXLnzmvoqhGZNYYmIiITIq1K0rokrUzC1bU1JkyYgbhx4xq6akRmj6GJiMhE7NmzEx06NFfjmCQkTZ06B40bNzN0tYgsBkMTEZGR8/X1xahRQzB16nh1PV++Ali8eC1cXLIZumpEFoWhiYjIiD169FCtvXTixFF1vW3bzhgzZrIa+E1EsYuhiYjISP31159qOYG3b98gQYIEmDVrIerWrW/oahFZLG7YS0RkhN1xgwb1RqNGdVRgKliwCA4dOsfARGRgbGkiIjIiDx7cV91xp04dV9c7deqOESPGw8HBwdBVI7J4DE1EREbUHdehQwu8e/dWrejt6bkYtWvXNXS1iOg/7J4jIjKi7jgJTIULF1PdcQxMRMaFLU1ERAaeHdeyZSOcPHlMXe/SpQc8PMbB3t6ezwuRkWFoIiIykF27dqB9e1c12DthwoTw8lrC1iUiI8buOSKiWObn54fhwwegQYOf/psdVxgHD55lYCIycmxpIiKL8fAh8Pp17D9u0qSAs/P/946T2XFHjx5S19u3d8Po0ZM4O47IBDA0EZHFBKZChQBv79h/bFkt4OxZ4Natv9GmTRO1d1z8+PExc+YC1KvXMPYrRERRwtBERBZBWpgMEZiEPO6kSXOxZEknaDQa5MmTD0uXrufecUQmhmOaiIhiweLFc1VgatGiLf755zgDE5EJYksTEVEscHBwxLRpS9C0aQuebyITxdBERBQqDQriHIriFHLjMpzwFV/hhMvIjVMoinMoCMBK53M3d+4y1KuXleeayIQxNBERBful6It2mA93TEM23EQArOAHW1hBAw2sYAs/WEODG3DBNLhjPtrBD3YRnsPMmRmYiEwdxzQREQm7AOTCRdWKNAtuyIpb//2S1MAevrCDn/pfrgv5vpST8rlxieeQyAIwNBERJf8GuzUH4eh1DI+LOqpWJW04CvuX5/cyEphOohjK4QDPI5GZY2giIouXOdE1+CYBzmfLitoTxqG4lxd2FS0aQWz6zg7+cIA3/kJNtjgRmTmGJiKCpY9hGoXB6usA6++/Es+6uKDGhAk6hycbBMAOvliOZup4RGSeGJqIyKLJoG8X3Ax2m7+NTaTDk4x5yo8L6nhEZJ4YmojIgmnULLmwRCU8dcd0dVwiMj8MTURksWQdJllWIKJB37qGJzlOdtxAAZyPwVoTkaEwNBGRxZKFK2UdJl3pEp7keHJcIjI/XNzSVPj5ARqN/EYmE+QTEMAOG0M/Bz4SjuyD3SYrfcvClZEVMjwVvXYNY+bPR5WzZ9Xx5LhEZH4YmkyBnx+sHz+EVYAvrL7/riYT4qPR4EaAL75xmItB3XrsqJakDEq2RpE1maJKG55O5ciBbl274kqrVmo8Uxx8iXZ9icj4MDSZAmlh8vMF7G0Aez5lpkZ2tv/mJ282K9hHoiuI9MvRLngrk5C95GRrlKiy8fdXwUnb0vSdFb4gTjRqSkTGip/ApsTGFrCLeI8rMjIa6VP1h72VNeytGJoMxc7ux193svmu7CUX1bBU6OZNjFy0CNVOnQqMXnI8OS4RmR+GJiKyWLJvXEQz53QNS1pyPDkuEZkfhiYisljnUBA34BLhDDpdwpKQ49yEC86jQIzVmYgMh6GJiCyYFabBHa3VgpRRD0tBTUd3dVwiMj8MTURk0eajHcpiR7TDki9scQm5MQ/tY7S+RGQ4DE1EZNH8YIfBGAXgHawDAtSmvZEJS8If1vCFHVyxAv78tUpktqIdmnx9fXHq1CkcPnwY9+/fx8uXL/H161ckS5YMyZMnR6FChVC2bFmkTZtWPzUmItKzO29zwO71ceR9eQsjFi1GrVPHde5gkxYmCUw18RcuIw+fGyIzFuXQtG/fPixYsACbN2/Gt2/fAtejCcnqvynWOXPmROvWrdG8eXMVqIiIjMYrR/j+Xg7ffJMiDTzV2k3y2yy8mXXawePSJSctTAxMROYv0qHpzz//xIABA3D16lUVkmxtbVGgQAEULVoUqVOnRpIkSeDk5IQ3b96oy5UrV1RLlPzfu3dvDBw4EO3bt8eQIUNUSxQRUWwI7Y+6YHytcQV5UQyn0A7z4Y5pajNfCUfft1qR+1updZgkTMksORn0LWOipIuPiMxfpEJTuXLlcOTIERWKGjZsiMaNG6N69epwdJTtCcJ3+/ZtrFmzBqtXr8asWbOwdOlSLFu2DL/88kt06k9EFKEvnz9j5tS5AGZEWFYCkBc6wwudUADn1ea7specbI0iK33LwpWyDtP3ZQU4S47IklhpIvzz6/+kFalbt25wd3dHokSJEJ2uvZEjR6JixYqqxckcfPjwAQkTJsT79++RIEECvR7b/5sv3v97FzZODrBx5F+0psZbE4ALvt6IwxXBDeLB/bvo1a0jbt+KK1vsGqYSAA4eBApw+SaiKPH3/75vfeLEwH9bPhrk8ztSLU0y0Dt+/PjRrZ8KS3L5+PFjtI9FRBSWg/v/weD+PfDp40ckSlQZ797xXBFR1FlHpnBogUmSWVTpI4AREYUUEBCAuZ7T4d6lrQpMBQoWwagJ03iiiCj2QlNoSpcujYcPH0b3MEREevHx4wf07Noec2d/D0kNf3fF3EUrkSlTQtjby+bJsc/BAUia1CAPTUTGtE6TzIorWbIktm3bpmbR6bKuk50dx+UQkf7duXUTvbp3wP17d2Fvb4+BQ0ejzq/11fdSp/HHpu2P8e5d5AdE+Pr64ZuvD7KmTa+OG1kSmJydI303IjK30DRo0CCMHj0a5cuXx9q1a1GjRo0wm8sXLlyIMWPG4O7du9F9WCKiYPb+vQtDB/TCly+fkSpVGkya7oVcefIFKyPBSS6R5ePjiy8+3sibWaNajYjIMkW7e05mwUkYkgUu69Spoxa8DEom5y1fvhzZs2dHx44d8eDBg+g+JBFRIH9/f3jOmIze3TuqwFSkaAmsXL/1h8BERGTw0CRatWqFHTt2IE6cOOjQoQMGDx6swtL69euRJ08etGzZUq3T5OLiotZmIiLSh48fPqCHWzssmDtLXW/i2hqeC5YjcRIOICIiI96wt3Llymr/udq1a2Ps2LGqxUn2oZPwlC1bNhWkmjRpAmtrveQ0IrJwMn6pZ7f2eHD/HhwcHDDYYyx++vlXQ1eLiMyY3kKTePv2LdKkSaO64F68eKEGfM+bNw+urq4MS0QUM+OXUqfB5BlzkTMXN8slopill2afAwcOoFKlSqhQoQKOHz+u9qCT1iWZKSeDw79+/aqPhyEiCycTSoKNXypWEivXbWVgIiLTCE0SlCQw7d+/H0mTJsXEiRPV+KUTJ06o7+3cuRNlypTBkydP9FNjIrLY9Zd+GL80fxnHLxGR6YSmgwcPqr1aRowYgTt37qBXr15qA1/Zx2XXrl1o2rQp/v33X5QoUQIXLlzQT62JyKLcvXMbzRvXxaEDe9X4pRFjJ6N3/yGwtdXrCAMiopgNTQMGDFDrLslA73jx4gX7noxpkuUGpMyjR49Qrlw5FaSIiHR1YN/fKjDJgpWy/tLC5etRu049nkAiMr3QJAtbJkqUKMIyMiD8y5cvai0nIiJdxi8tmDNTbYny+fMnFCpSDCvWbUGu3Hl58ojIIGJt/n/btm2xdetW1bRORBSeL58/o2+PLvCcOUUtWyL7x3ktWIEkSZPxxBGRwcTqgADZYuXQoUOx+ZBEZGIePriPXt064NbN67Czs8eAISNQ97dGhq4WEVHstTRp5c+fP1Yex9PTE5kyZVKD0gsXLhxhWJNlE6SclM+cOTPmzJkTK/Ukov87cewwXBv9ogJTsmTJMW/JagYmIjLNlqbOnTtj4MCBSJcuXbQfeM2aNWrMgqwSrm+yNpS7u7sKTqVLl8bcuXNRs2ZNXLlyBenTp/+hvAxkr1WrFtq1a4cVK1bgyJEj6mdNnjw5fvvtNxiDz1+/wEbjD5sAux++Z2NjDUcHx/+X/fIlzONYW1vBydEpSmW/fP2qukpCY2VlhThOUSv79dtXBASEXlbEjRMnSmW/eX+Dv3+AXspKfaXewtvHG35+/jqX/ezri6++3rCysob/f7drOTg6Bi786uvjAz8/vzCPG5my9g4OsLGxiXxZX1/4+fqGWdbO3j5wxlpkysrjSz3CYmNri3Wrl2PapDHq90Ku3HkwZsIMJE+REl9DvEZt7ezUJBPtvnM+3t5hHlceX+oR2bJSB+9v3wK/5yPPoTyXnz+rn0Ue3z5I2fDWopPjaoclyHtCxnbqo6w8Z/JHnpbUTR9l5TXmFOT9GZmyUt9w3/dB3nORKSvnV85zWOLGjRulsrJnqrwu9FFW6hv4vvf2Dvc9F5mycn6173sfHx/1vtNHWXk9aN/3kSnr6+uryodFXr9Bf0foWlbOgZyLsNjY2MPK6sfPv1iniQRra2uNg4ODpk2bNpoDBw5oIuvFixeamTNnanLkyKGONWrUKE1MKFasmKZjx47BbpPH7N+/f6jl+/btq74fVIcOHTQlSpQI8zG+ffumef/+feDl4cOH8u5XX+ub31cfdeywLtXKV9J8vP448BLHySnMsmWKlQxWNmniJGGWLZQnf7Cy6dOmC7NsjqzZgpWV62GVleMELSuPE1ZZqV/QslL/sMrKzx20rJyX8M5b0LJ1q/8Ubtln524Glm3ya4Nwy945diGwbLsmLcItu233Ic3Zy3fVxbVVu3DLrt+yK7Bs+87dwy27fM3mwLLde/UPt+y8xasDy/Yb5BFu2emeCwPLDh81Mdyy46fMCiwrX4dXNn+BwoFfFytROtyyUkftcaXu4ZWVn11bVs5JeGXlnGrLyrkOr2y3br01Hz5o1OXixbvhlm3btnNg2Tt3XoRbtkmTFoFlnz79FG7ZunXrB5aVS3hlq1WrFaxsnDhxwixbpkz5YGWTJk0WZtmCBYsEK5s+fYYwy+bIkStYWbkeVlk5TtCy8jhhlZX6BS0r9Q+rrPzcQcvKeQnvvAUtK+c7vLLyfGnLyvMYXll5HWjLyusjvLLy+tKWldddeGVPnLgUWLZ//2Hhlt2372Rg2ZEjJ4Rbdvv2fYFlJ00K/728bt22wLJeXovDLbt06brAsvJ1eGVnz16seflSo/Hz0/tHrPrc1vXzO1ItTbLOUt++fbFo0SIsXrxYbZkiLTjFihVTXVuyEniSJEnUX2Dv37/HmzdvcPXqVZw+fVrtSycLYEpSl0Uwp06dik6dOuk9BEqqPXPmDPr37x/s9mrVquHo0aOh3ufYsWPq+0FVr14dCxcuVElZ+1dtULK/noeHh55rT2SZ/j1/Rv0l26PPIKRzTo+Tx48YukpERD+wkuQEHUnXnKy1JOsxzZw5E1u2bFEhRdvMGBbtQ2TNmlV1gXXs2BHx48dHTJCVx9OmTau62EqVKhV4+5gxY7B06VJcv379h/vIli8tW7ZUP5+WBCzp2pPjSRgMSZoRgzYlfvjwAc7OziosymKf+uT/zRdPjl+GjaMDbBzZPWeK3XOXfL0Rx8oa9uyeC3TxwnkM6tsdr1+/QoIECTF+6mwUL1E6wq48Q3XPffHxRu6MWVWXArvnvmP3XJD3PbvnYqV7LnFi+Rp6JZ/fsiC3Lp/fkWppGjdunAojMmhatkZ59eqVanGSFigJGQ8fPgzWNysPXqBAAVVWWnIkcMWWkEFOglt44S608qHdHvTJjs3lE+I6xYGNU+ihKbwwoM+yQcch6bNs0HFT+iwbdJyXPss62DvAwV73srCzg5OvDZxCCU1ByQe29kM7IjFWNkgg0WdZ+cUYdPXuLZvWY4zHYPj6+iCrS3ZMmTlPtTCFVjY80jrlpONrODJlVRgIUtbGxxcaGxs1xiXk+17KBh37Eh75fRITZYUxlA06DkmfZYOOm9Jn2aDjvPRZNjKfD5EpK7042rF0hiprF0u/I0KS4WThDP2KNZEKTfIDBR0IlyxZMvTr1w+1a9dW+82Jd+/eqQFz2m662CZ1kl+Oz549C3b7ixcvkDJlylDvkypVqlDLy88rXYlEpB/yR9XUiaOxesUSdb1ileoYOWYy4kTig5mIyCSWHEiRIoWaaRayR0+CkpasDi4hxBCBScjjyviqPXv2BLtdrgftrguqZMmSP5TfvXs3ihQponNKJqLwvXv3Fm4dWgYGpg5d3DFxqicDExGZZ2gqX768apFp0aIFzp8/H+Z0UUPr2bMnFixYoAasy0D0Hj164MGDB2oslZC98Jo3bx5YXm6/f/++up+Ul/vJIPDevXsb8KcgMi4+Gk2U3/O3b91Q+8fJAG8npziYOM0LHTp3D5weTURkCiLVPTdy5Ej8888/ai2jlStXBvZJS+uTDAzPly+fuiSWkVoG1KhRI7x+/RojRozA06dPkSdPHuzYsQMZMmRQ35fbJERpySKY8n0JV7Nnz1azAmfMmGE0azQRGdozfz80e/MEqWxs0TluYpS0d4xwAojWgb17MKhfD3z58hlp0zmr8Usu2XLEeJ2JiAw6e068fPkS06ZNU/vIyWKRoQ2ylsUvZeXvoBcXFxeYs8iMvo/K7Ln3/97VeSA4GRdvTQAuhDF7zlRc9fVG07dPIbWXXxi5bO0jDE/ye2HRfE94zpisvi5StISaIZc4cRKYGh+f77Pn8mZ24f6ZRAagHQhu6NlzkQ5NQX38+FE9kGw7UqdOHfz7779qJp208gQ+wH+/UGVGQ968eVWAMsctShiayBJCk5Z0qgWEE55kVWaPIX2x+69t6rpsuNur3xCTHSPI0ERkWMYSmqK1Ya92rSXpzpoyZUrg7Y8fP1YBKujl5s2bOHHiBE6ePGmWoYnIkmhXtLrm5wO398+DhacXz5+hR9f2uHblkpqB2m+QB35rqP/tkoiIYlu0QpO4c+dOsJYlIYtLykX2c9OSZQguXryoAhQRmWd4yuTjh1fjx+HjlUtIlDgJJk7zROEixQ1cSyIiIwlNGTNmVBddFgYrWrSouhCReYanuzZWwJDBcHBthh7pMqNQamcD14yISH8435eI9Ee7E3rGjBhmGwDXt09x1Pur0S5PQkQUGQxNRBTj3XYSnk74fOWZJiKTxtBERDEenq74+WDCxzc800Rk2WOaiIjCEnRpgq7xDLvoLRFRdDE0EVGMhaUcOiyCSURkKhiaiEhvGJaIyJwxNBGRTqv/h4dhiYgsAQeCE1G47ty6iaEDeob7C0S64WYlTInliVOjlIMTu+KIyCyxpYmIwnT44D4M6N0Nn9OkDnY7W5aIyBIxNBHRD2QxypVLF2LqpDHq6xy5cuOabMAt3+MAbyKyUAxNRBSMr48Pxowcgi0b16nrdX9riNYDh6PVx5dIaWPL2XBEZLEYmogo0Ns3r9HbvRPOnTkFa2tr9OgzCE1cW6kxStsdnGEnrU1cOoCILBRDExEpt2/dQPfObfDk8SPEixcfYyfNQOmyFQLPjj3DEhFZOM6eIyIcOrAXLZv8pgJTOucMWLJqY7DAREREDE1EFk0Gea9YsgDuXdri8+dPKFK0BJat2YTMWbIaumpEREaH3XNEFirkgO9f6zdG/0EesLO3N3TViIiMEkMTkQV6+/YN+rh3wtnTJ38Y8E1ERKFjaCKywBW+3d3a4tHDB6EO+CYiotAxNBFZkKOHD6B/r6749Okj0qZzxrTZC5AlazZDV4uIyCRw9hyRhVizcim6dWqtAlOBQkWwbPUmBiYiokhgSxORmfP19cWkcSOwfs0Kdf3nX37DoOGjYW/vYOiqERGZFIYmIjP28cMH9OvVBcePHlaDvLv26IsWrTtwwDcRURQwNBGZqYcP7qN7lza4d+c2HJ2cMHr8NFSsXM3Q1SIiMlkMTURm6MzpE+jTvRPevXuLlKlSY+qs+ciRM7ehq0VEZNI4EJzIzGzdtAGd2riqwJQ7Tz4sW7OZgYmISA/Y0kRkJgICAjBr2kQsWThHXa9avRY8xkyGo6OjoatGRGQWGJqIzMDXL18wuH8P7Ptnt7rermNXdOjirlb7JiIi/WBoIjJxL54/UxvuXrt6GXZ29hg2ajxq1a5r6GoREZkdhiYiE3bl8kX0cGuHly+eI3GSpJgyYy7yFyxs6GoREZkltt0Tmai9f+9C2xaNVGDKnMVFrfDNwEREFHMYmohMjEajUYO9+7h3wrevX1GydDksXrlB7SVHREQxh91zRCbE18cHY0YOwZaN69T1hr+7onf/obC15VuZiCim8TctkYl4/+4d+vTojNMnj6lZcb37D0Hjpi0NXS0iIovB0ERkAh7cv4vundvg/r27iBs3HsZOmoEy5SoaulpERBaFoYnIBLZE6d2tI96/f4dUqdNguudCuGTLYehqERFZHA4EJzJif27+Q22JIoEpT978WLZ6MwMTEZGBsKWJyEi3RPGaOQUL581W16vW+AkeoydxSxQiIgNiaCIyMt++fcOwgb2wZ9cOdb1N+y7o1LUnt0QhIjIwhiaimKDRwPrKJVhfPA+bm9eR/etnWDs4wcclG77kyYevOXMDVlY/3O31q5fo0bU9Ll04D1tbOwzxGIuf6/7G54iIyAgwNBHpk68vbNethP3SBbC+fxcaCUY2NkilAWAFWPn7w0qjwbcMGfHStTVe1W8M2Nmpu966eV3NkHv65DESJkyESdO9ULhoCT4/RERGgqGJSA+8NQFwvHkDjr27wvr6lcBWJAlI8POTvBSMw4P7SDdqKJKuX43746di7/Nn6N+rKz59+oj0GTJihtcipM+Qic8NEZERYWgiiqZH/r4o/+Q20j++iZGJ4qO6NiyFQ/t9p1s3kKX+z+gS4I9PAQEoVKQYJk2fg0SJEvN5ISIyMlxygCiaXt+9jRc21jjr4oKa48ahuJcXdhUtivBj03fSXWfj54ttAQHoXKEyPOcvY2AiIjJSDE1E0eHrC/upE9SXAdbf304SnmpMmKBzeLKR7jorK0x8+gT2VnxLEhEZK/6GJooGGfQtA76D8rexiXR4stFo4HT9KpJtWMPng4jISDE0EUWVRqNmyf0wyjuq4cnKCslXLFbHJSIi48PQRBTVN4+sw3T/boSDvnUNT3Icx3t34XTtCp8TIiIjxNBEFNU3z8Xz39dh0pEu4UmOF+fiv3xOiIiMEJccIIoi61s3AJvIv4VChqei165hzPz5qHL2LDQ2NrC/dRM+7KIzKj46zYUkInPH0EQURVbfvknbUJTPnzY8ncqRA926dsWVVq3U4QK+fcEXTQCfF6OigaOs6G7oahCRQTE0EUWRxtHx+94oUWTj76+Ck7alSUhvX+I48ZDPzoHPi1GxgZU/YPffshJEZJkYmoiiKCBrNsDfL8phqdDNmxi5aBGqnTr1/+gle9O5ZIcD12syLpJmrazYSUdk4RiaiKIoIG+BCGfO6RyW/iPHk+MSEZHxYWgiiqKAXHkQkCFThDPodAlLQo6jyZAJATlz8zkhIjJCDE1EUWVlBZ8WbYHVS6IVlgJpAJ/mbb53BRERkdFhaCKKho916uHJlnXRC0uSl2xsEOCSHX6NmvH5ICIyUgxNRFH08vUrNOzUEm++fFHXrQMC1Ka9kQlLQiMzsuzs4D1xJmDLtyQRkbHib2iiKLh2+yYatG+Oe48eIGHmLEji7YPM9+5i5OLFqH7ihM4LEUgLkwSmr/NXICBbDj4XRERGjKGJKJIOnTiKJm5t8e7De2RyzoANsxchQ8bMcPS1heP7j2rtJo3MUA9nZp0aPC4LWbpkVy1MDExERMaPK7URRcLqzRvwS5smKjAVK1gY/6z7E9kyZ1XrKmmy5cDXDTvgPXQUNOkz/n9GnK0dNLa23///b5C3zJKTclKegYmIyDSwpYlIBxqNBmNnTVEX8WuN2pg7fhqcHJ2CF7Szg1/TlvBr0gLWVy+rTX2tb96A1bev0Dg6IcAlm1qHSS0rwFlyREQmhaGJKAI+Pj5wG9wbq7f8oa73aNcFw3v2h3V4W2pYWX1fxylXHp5fIiIzwdBEFI6379+hqVtbHDp5DDY2NpgybAxac1kAIiKLxNBEFIZ7Dx/gt/auuHHnFuLHjYdl0+eiStkKPF9ERBbK7AaCv337Fq6urkiYMKG6yNfv3r0Ls7yvry/69euHvHnzIm7cuEiTJg2aN2+OJ0+exGq9ybicvnAOlRr9rAJTmpSpsGvVRgYmIiILZ3ahqUmTJjh//jx27typLvK1BKewfPnyBWfPnsWQIUPU/xs3bsSNGzdQp06dWK03GY8//96JWq711eKV+XLmxr7125A3B/eDIyKydGbVPXf16lUVlI4fP47ixYur2+bPn4+SJUvi+vXryJ49+w/3kdaoPXv2BLtt5syZKFasGB48eID06dPHWv3J8GYvmY8B4zzUbLlq5SphyVQvxI8Xz9DVIiIiI2BWoenYsWMqBGkDkyhRooS67ejRo6GGptC8f/8eVlZWSJQoUZhlvL291UXrw4cP0aw9GZK/vz8GjPWA1/KF6nqbxq6YNGQUbLmtCRERmWP33LNnz5AiRYofbpfb5Hu6+PbtG/r376+6+RIkSBBmubFjxwaOm5KLs7NztOpOhvP5yxe1wrc2MI3oPQhTh49lYCIiItMLTcOHD1ctP+FdTp8+rcrK1yFJV0tot4c2KLxx48YICAiAp6dnuGUHDBigWqS0l4cPH0bjJyRDefHqJX5q3gA79u6Gg70Dlk2bgx7tOuv0eiEiIstiEt1zbm5uKsyEJ2PGjLhw4QKeP3/+w/devnyJlClTRhiYGjZsiLt372Lv3r3htjIJBwcHdSHT3nS3fjtX3H/8EEkSJcZar8UoUaiooatFRERGyiRCU7JkydQlIjLgW1p9Tp48qQZyixMnTqjbSpUqFWFgunnzJvbt24ekSZPqtf5k3JvuZsmQEX/MX4EsGTIZulpERGTETKJ7Tlc5c+ZEjRo10K5dOzWDTi7yde3atYMNAs+RIwc2bdqkvvbz80P9+vVV997KlSvVgGAZ/yQX2T6DzM+aLX8EbrpbvGAR/L3mTwYmIiKyrNAkJPjIQpXVqlVTl3z58mH58uXBysjyA9L6JB49eoStW7eq/wsUKIDUqVMHXmTGHZkPGds2wXMa2vXtploXZdPdbUvXIlmSJIauGhERmQCT6J6LjCRJkmDFihURfngGHQsV9DqZJwlJ3Yf1x/I/1qjr3dt0wojeA8PfdJeIiMicQxNRSB8+fYRrt/bYe+SgCkmTBo9Eu6YteaKIiChSGJrIrD1+9gT12zfHpetXEcfJSa3wXbNiVUNXi4iITBBDE5mti9cuq8D05PkzpEyeAuvnLEXBPPkMXS0iIjJRHNBBZunvQ/tRvUk9FZhyZM2GvWv/ZGAiIqJoYWgis7N0/WrU79AcHz9/QtliJbFn9WakT5vO0NUiIiITx9BEZkNmQY6cNgFug3ur9bYa16mHTQtXIlGChIauGhERmQGOaSKzIAuRSlhaveUPdb1vp+4Y3L0P95AjIiK9YWgikycrezd1a4uDJ47CxsYG0z3Go0WD3w1dLSIiMjMMTWTSHj55jN/au+LqzeuIFyculs+YhyplKxi6WkREZIYYmshk/XvlklpS4NnL50idIhU2zFuKfDnzGLpaRERkpjgQnEzS7gN7UaNpPRWYcmXLgb3rtjIwERFRjGJoIpOzZP0qNOzUEp++fEaFkmWwe9UmpEud1tDVIiIiM8fQRCa1pMCIqePRdXAftaTA73Xr4495y5EwfgJDV42IiCwAxzSRySwp0GVQL6zZulFd79fZHYO69eaSAkREFGsYmsjklhSYMWICmtdvbOhqERGRhWFoIqPGJQWIiMhYMDSR0bpw9RJ+a8clBYiIyDhwIDgZpb8P7Uf1Jt+XFMjpkh3/rOWSAkREZFgMTWR0lm1Yg/odmqslBcoVL6WWFHBOwyUFiIjIsBiayKiWFBg9Y5KaJSdLCjSuUw+bFqxEogQJDV01IiIijmki4+Dr64uuQ/ti5cZ16nqfjt0wxL0vlxQgIiKjwYHgZHAfPn2Ea7f22HvkoFpSYMqwMWjdqJmhq0VERBQMQxMZ1JPnT9WmuxevXUHcOHGwdNocVC9fmc8KEREZHYYmMpirN6+jXrtmePT0CVIkS44Nc5ehYJ58fEaIiMgocSA4GcTB40dQ9fe6KjBly5wV/6zZysBERERGjaGJYt3arRtRt21TvP/4ASULF8Oe1ZuR0Tk9nwkiIjJqDE0Uq0sKTJo7E237dFWz5erV/BlbF69GkkSJ+SwQEZHR45gmihV+fn7oNWIQFq1doa53bd0Bo/oMhrU1czsREZkGhiaKcZ+/fEHLnp2wc9/fat2l8QM90Kl5G555IiIyKQxNFKNevHqJBh1a4Oylf+Ho4IgFk2bil2q1eNaJiMjkMDRRjLl59zbqtW2Ge48eqHFLa70Wo0ShojzjRERkkjighGLE8bOnUKXxLyowZXLOgL/XbGFgIiIik8bQRHq3ZfcO/NyyMd68e4vCeQvgn7Vb4ZIpC880ERGZNIYm0iuvZQvVPnLfvL+hRsUq2L5sPZInTcazTEREJo9jmkgvAgICMHjiKMxcNFddb9PYFZOGjIKtLV9iRERkHviJRtEmrUod+rlj419/quvDew1Az3Zd1PICRERE5oKhiaLl7ft3+L1LGxw5dRx2dnbwHD0ZjX/5jWeViIjMDkMTRdmDx49Qr10zXL99EwnixcfKWfNRoWRZnlEiIjJLDE0UJReuXsJv7Zrj2cvnSJMyFf6Ytxx5cuTi2SQiIrPF2XMUaf8cPoDqTeqpwJQrWw7sXfcnAxMREZk9hiaKlBUb16J+h+b49OUzyhUvhV0rNyJtqjQ8i0REZPYYmkgnGo0GEzynodOAnvDz80PD2r9i44IVSJQgIc8gERFZBI5poghJSOrhMRBL1q1U13u2d8OwHv1gbc3MTUREloOhicL1+csXNHfvgN0H9qqQNGnwSLRr2pJnjYiILA5DE4XpxauXaNChBc5e+hdOjo5YNNkTtatU5xkjIiKLxNBEobp59zbqtW2Ge48eIGniJFg3ZwmKFSjMs0VERBaLg1LoByfOnUaVxr+owJQ5fUb8vWYLAxMREVk8hiYK5s+/d6J2i0Z48+4tCuctoAJT1oyZeZaIiMjiMTRRoHkrl6CpW1u1AW+NilWwfdl6JE+ajGeIiIiI3XMkAgICMHTSGPQaMUitx9S6UTOsnrUQcePE4QkiIiL6DweCWzhvH290HtAL67ZtUteH9uiH3h26wsrKytBVIyIiMioMTRbs/ccPaNKlDQ6eOApbW1vMHj0JTeo2MHS1iIiIjBJDk4V69PQxfmvfHFduXEP8uPGwYuZ8VCpdztDVIiIiMloMTRbo8vWrqNeuGZ48f4ZUyVPij/nLkC9nHkNXi4iIyKhx9pyFOXDsMKo1+VUFpuxZXLB33Z8MTERERDpgaLIg6/7chF/bNcOHTx9Rqkhx7Fm9Gc5p0hq6WkRERCaBockCyDIC0xZ4oU1vN/j6+uLXGrWxZdEqJE6YyNBVIyIiMhkc0xSLHjwAXr2K/P38vYFPNxxg7WgPG/vIPWX+Af6Yt3ISVm2arK53adEWY/oPg7U18zIREVFkMDTFYmDKnh349i0q97YDkCEajz4SwHqMHdAKbi3bR+M4RERElouhKZZIC1PUApM+OGFErylwa1nSUBUgIiIyeeyjsRAVuQYTERFRtDA0EREREemA3XNGS4OCOIeiOIXcuAwnfMVXOOEycuMUiuIcCgLg/nBERESxhaHJyNjCF+0wH+6Yhmy4iQBYwQ+2sIIGGljBFn6whgY34IJpcMd8tIOfGihOREREMYndc0YkNy6pVqRZcENW3FK3SUCyhy/s4Kf+l+tCvi/lpLzcj4iIiGIWQ5ORKIcDOIliKgBJMNKGo7Boy0h5uZ/cn4iIiGIOQ5MRkODzF2rCAd6wg3+k7ivl5X5yf7Y4ERERxRyGJiMYw7QCzWAHH9ggIErHkPvZwRfL0Uwdj4iIiPSPocnAZNB3PlyIdAtTSDLmKT8uqOMRERGR/pldaHr79i1cXV2RMGFCdZGv3717p/P9O3ToACsrK0ybNg0xT6NmyelTd0xXxyUiIiL9MrvQ1KRJE5w/fx47d+5UF/lagpMuNm/ejBMnTiBNmjSIDbIOkywrENGgb13JcbLjBgrgvF6OR0RERGYamq5evaqC0oIFC1CyZEl1mT9/PrZt24br16+He9/Hjx/Dzc0NK1euhJ1d7Kx7JAtXyjpM+iTHk+MSERGRfpnV4pbHjh1TXXLFixcPvK1EiRLqtqNHjyJ79uyh3i8gIEC1RvXp0we5c+fW6bG8vb3VRevDhw+Rrq+s9C0LV8r6S/oix5PjhuTv4wf/bxwkThQl/n6QYYcaf3P7U5PINAREbZ6U3plVaHr27BlSpEjxw+1ym3wvLOPHj4etrS26deum82ONHTsWHh4eiA7ZGkVW+tYvDeLgyw+3Bnzzgf/X/4c8IookWzsE+Ftx9yIiA7G1BawMvHuYSYSm4cOHRxhQTp363iUlg7hD0mg0od4uzpw5g+nTp+Ps2bNhlgnNgAED0LNnz2AtTc7OzogM2UtOtkbRLyt8QZwfbo2XLQ0S5tfzQxFZEvn9IL+1ichgb0FrA7f0msRvABlr1Lhx43DLZMyYERcuXMDz589/+N7Lly+RMmXKUO936NAhvHjxAunTpw+8zd/fH7169VIz6O7duxfq/RwcHNQlOmTzXdlLTp/keHLckGwc7GDjqNeHIiIisigmEZqSJUumLhGRgd/v37/HyZMnUaxYMXWbzIaT20qVKhXqfWQsU5UqVYLdVr16dXV7q1atEJNk3zh9zZzTkuPJcYmIiMgCQ5OucubMiRo1aqBdu3aYO3euuq19+/aoXbt2sEHgOXLkUGOSfv31VyRNmlRdgpLZc6lSpQpz4Li+nENB3ICL2nxXH+FJZs7dhAvOo4Be6kdERET/Z3bzQGTJgLx586JatWrqki9fPixfvjxYGVl+QFqfDM8K0+Cu1yNOR3d1XCIiItIvK42MkqZok4HgsrSBhLEECRL88P2zZ4HChX+8n+wVJ91pstludLZS8YUtLiE3iuI0/ENpQDxzBihUKMqHJyIissjPb7NuaTI1frBDM6yAL+zhH8WnQ+7nCzu4YkWogYmIiIiij6HJCFxGHtTEX/CGg2oxigwpL/eT+8txiIiIKGYwNBmJgyiPYjiputhkQHdE26toy0h5uZ/cn4iIiGIOQ5MRkZaiYjgFN8zCLWRVt0kw8oEdfGCr/teGKZklJ+WkPFuYiIiIYh4HwBjhGCcvdIYXOqEAzqvNd2UvOdkaRVb6loUrZeD492UFOEuOiIgotjA0GS0rnEdBdSEiIiLDY/ccERERkQ4YmoiIiIh0wNAUS2TrPEcDbZgrj6vD1n1EREQUDo5p0hPtwuqysmhoEiUCTp8GXr9GrJOt9eTxw6gaERGRxfrw34ejLhukMDTpycePH9X/zs7O+jokERERxeLnuGynEh7uPacnAQEBePLkCeLHjw8rKyu9p2AJYw8fPoxwXxzieTZ2fD3zPJsTvp5N/zxLC5MEpjRp0sDaOvxRS2xp0hM50enSpUNMkhcKQ1PM43mOHTzPPM/mhK9n0z7PEbUwaXEgOBEREZEOGJqIiIiIdMDQZAIcHBwwbNgw9T/xPJs6vp55ns0JX8+WdZ45EJyIiIhIB2xpIiIiItIBQxMRERGRDhiaiIiIiHTA0ERERESkA4YmI+Dp6YlMmTLB0dERhQsXxqFDh8Itf+DAAVVOymfOnBlz5syJtbpa0rneuHEjqlatiuTJk6vF1EqWLIldu3bFan0t5TWtdeTIEdja2qJAgQIxXkdLPM/e3t4YNGgQMmTIoGYhZcmSBYsWLYq1+lrKeV65ciXy58+POHHiIHXq1GjVqhVeG2LjURNy8OBB/Pzzz2pVbtlVY/PmzRHexyCfhRoyqDVr1mjs7Ow08+fP11y5ckXTvXt3Tdy4cTX3798PtfydO3c0ceLEUeWkvNxP7r9hw4ZYr7u5n2v5/vjx4zUnT57U3LhxQzNgwAB1/7Nnz8Z63c35PGu9e/dOkzlzZk21atU0+fPnj7X6WtJ5rlOnjqZ48eKaPXv2aO7evas5ceKE5siRI7Fab3M/z4cOHdJYW1trpk+frn5fy/XcuXNr6tatG+t1NyU7duzQDBo0SPPHH3/IrrmaTZs2hVveUJ+FDE0GVqxYMU3Hjh2D3ZYjRw5N//79Qy3ft29f9f2gOnTooClRokSM1tMSz3VocuXKpfHw8IiB2pmPqJ7nRo0aaQYPHqwZNmwYQ1MMnOe//vpLkzBhQs3r1691OTxF8TxPnDhRhf+gZsyYoUmXLh3PqY50CU2G+ixk95wB+fj44MyZM6hWrVqw2+X60aNHQ73PsWPHfihfvXp1nD59Gr6+vjFaX0s716FtyiybOiZJkiSGamm553nx4sW4ffu2WryOYuY8b926FUWKFMGECROQNm1aZMuWDb1798bXr195yvV4nkuVKoVHjx5hx44daiPY58+fY8OGDfjpp594nvXIUJ+F3LDXgF69egV/f3+kTJky2O1y/dmzZ6HeR24Prbyfn586nvSfk37OdUiTJ0/G58+f0bBhQ55iPb6mb968if79+6txIjKeiWLmPN+5cweHDx9W4z82bdqkjtG5c2e8efOG45r0eJ4lNMmYpkaNGuHbt2/qd3OdOnUwc+ZMvrT1yFCfhWxpMgIy6C0o+esk5G0RlQ/tdor+udZavXo1hg8fjrVr1yJFihQ8tXo6z/KB1KRJE3h4eKiWD4q517O0lMr35AO9WLFiqFWrFqZMmYIlS5awtUmP5/nKlSvo1q0bhg4dqlqpdu7cibt376Jjx46Re3IpQob4LOSfdQaULFky2NjY/PAXy4sXL35I0FqpUqUKtbz8hZ40adIYra+lnWstCUpt2rTB+vXrUaVKlRiuqWWdZ+nulOb0c+fOwc3NLfDDXX75yWt69+7dqFSpUqzV35xfz/KXt3TLJUyYMPC2nDlzqnMt3UkuLi4xXm9LOM9jx45F6dKl0adPH3U9X758iBs3LsqWLYtRo0axN0BPDPVZyJYmA7K3t1fTJffs2RPsdrkuTbyhkWnvIcvLB4uMVbCzs4vR+lrauda2MLVs2RKrVq3imIQYOM+ylMPFixdx/vz5wIv8RZ49e3b1dfHixXV/ki1IVF7P8kH+5MkTfPr0KfC2GzduwNraGunSpYvxOlvKef7y5Ys6p0FJ8AraEkLRZ7DPwhgdZk46T2dduHChmjbp7u6uprPeu3dPfV9maLi6uv4wzbJHjx6qvNyPSw7EzLletWqVxtbWVjN79mzN06dPAy8yNZ70d55D4uy5mHk9f/z4Uc3gql+/vuby5cuaAwcOaFxcXDRt27bly1mP53nx4sXq94anp6fm9u3bmsOHD2uKFCmiZuFR2OT1ee7cOXWRaDJlyhT1tXZpB2P5LGRoMgLyoZwhQwaNvb29plChQuqXmVaLFi005cuXD1Z+//79moIFC6ryGTNm1Hh5eRmg1uZ/ruVrefOGvEg50t95DomhKWZez+Lq1auaKlWqaJycnFSA6tmzp+bLly98Oev5PMsSA7I8iZzn1KlTa5o2bap59OgRz3M49u3bF+7vW2P5LLSSf2KuHYuIiIjIPHBMExEREZEOGJqIiIiIdMDQRERERKQDhiYiIiIiHTA0EREREemAoYmIiIhIBwxNRERERDpgaCIiIiLSAUMTERERkQ4YmoiIiIh0wNBEREREpAOGJiIiIiIdMDQRERER6YChiYgohIEDB8LKygrlypUL9dwMGTJEfT9Xrlx4+/Ytzx+RhbDSaDQaQ1eCiMiYSBDKlCkT3r9/j71796JixYqB31u8eDFat26NVKlS4dixY8iYMaNB60pEsYctTUREISROnBi9evVSX3t4eATe/s8//6BDhw6IGzcutm3bxsBEZGHY0kREFIqPHz8ic+bMePXqFfbv34+kSZOiTJky+PTpE7Zs2YKffvqJ543IwrCliYgoFPHjx0e/fv3U1/3791chSbrrZs+ezcBEZKHY0kREFIavX78iS5YsePr0aWB4Gjt2LM8XkYViaCIiCsO7d+9QvHhx3LhxA8mTJ1fhycbGhueLyEKxe46IKBQ+Pj6oV6+eCky2trZ4+fIlNmzYwHNFZMEYmoiIQtGmTRvs27cP1apVw7x589Rtw4cPh7+/P88XkYViaCIiCmHQoEFYsWIF8ubNi/Xr16N58+bIli0brl27pm4nIsvEMU1EREEsWLAA7dq1Q5o0aXD8+HE4Ozur25cvX67CkyxDcP36ddVlR0SWhS1NRET/2blzJzp16oR48eKpxSu1gUk0adJEtTbduXMHixYt4jkjskAMTUREAM6fP48GDRpAdpZau3YtChYsGOy8yKw52XNOjBo1Ct7e3jxvRBaG3XNEREREOmBLExEREZEOGJqIiIiIdMDQRERERKQDhiYiIiIiHTA0EREREemAoYmIiIhIBwxNRERERDpgaCIiIiLSAUMTERERkQ4YmoiIiIh0wNBEREREpAOGJiIiIiIdMDQRERERIWL/A3z22n1ZTJMKAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import nonlinear_plots\n",
    "nonlinear_plots.bisection_root()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr />\n",
    "\n",
    "**Exercise**: Apply the bisection algorithm *by hand* to solve $f(x) = x^3 + 2x^2 - 5x - 6 = 0$ for $x\\in [1, 4]$.  Find the root to within $\\tau = 10^{-1}$.  Write out all values of $L$, $R$, $C$, and $f(C)$ for each step.\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**: Implement the bisection algorithm in Python  and use it to solve $\\sin(x) - 2/5$ for $x \\in [0, 1]$.\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**: Modify the bisection algorithm to accept a desired number of iterations $n$ instead of the tolerance $\\tau$.  Implement this modified algorithm in Python.  Then, show the error in the approximate root of $f(x) = \\sin(x) - 2/5$ for $x\\in [0, 1]$ as a function of $n$.\n",
    "\n",
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Newton's Method\n",
    "\n",
    "A more powerful way to find roots of $f(x) = 0$ is [Newton's method](http://mathworld.wolfram.com/NewtonsMethod.html), sometimes called the Newton-Raphson method.  Like bisection,  Newton's method produces a sequence of approximations for a root.  The values of the sequence are increasingly close to the root.  Unlike bisection, Newton's method requires not a range in which a single root lives but an initial guess for what the root is.  \n",
    "\n",
    "Let the root of interest be $x_r$.  Assume we have a good way to guess the root (intuition, a graph, etc.), and call that initial approximation $x_0$.  Then, for some $\\Delta$, $x_r = x_0 + \\Delta$,Taylor series provides the following relationship:\n",
    "\n",
    "$$\n",
    "   f(x_0) + \\Delta \\frac{df}{dx}\\Bigg |_{x=x_0} + \\mathcal{O}(\\Delta^2) = 0 \\, .\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formally, this equation is still exact because $\\mathcal{O}(\\Delta^2)$ captures *all* terms proportional to $\\Delta^2$ and higher powers of $\\Delta$.  However, if $\\Delta$ were small enough that the $\\Delta^2$ terms can be ignored (i.e., our guess $x_0$ were close enough to $x_r$), then we're left with the *approximate* relationship\n",
    "\n",
    "$$\n",
    "  f(x_r) = 0 \\approx f(x_0) + \\Delta f'(x_0) \\, ,\n",
    "$$\n",
    "\n",
    "from which it follows that\n",
    "\n",
    "$$\n",
    " \\Delta =  -\\frac{f(x_0)}{f'(x_0)} \\, .\n",
    "$$\n",
    "\n",
    "We expect $x_1 = x_0 + \\Delta = x_0 - f(x_0)/f'(x_0)$ to be a *better* approximation than $x_0$.  Moreover, the same process can be applied to determine a sequence of estimated roots $x_2$, $x_3$, and so on.  This process is Newton's method, and is presented in the following pseudocode:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "'''Newton's method for finding the root f given x_0'''\n",
    "Input: f, fp, x_0, tau\n",
    "Set x = x_0\n",
    "While |f(x)| > tau do\n",
    "   # Compute the Newton \"step\" \n",
    "   Set Delta = -f(x)/fp(x)\n",
    "   # Compute the next value of x\n",
    "   Set x = x + Delta\n",
    "Output: x\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr />\n",
    "\n",
    "**Exercise**:  Consider $f(x) = \\cos(x) = 0$.  With $x_0 = 1$, apply *one* step of Newton's method.  What is the result (i.e., $x_1$)?\n",
    "\n",
    "*Solution*: $x_1 = x_0 - \\frac{\\cos(x_0)}{-\\sin(x_0)} \\approx 1.642$\n",
    "\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**: Apply Newton's method *manually* to solve $f(x) = x^3 + 2x^2 - 5x - 6 = 0$ to within $\\tau = 10^{-5}$. Use $x_0 = 3/2$. Write out $n$, $x_n$, $f(x_n)$, $f'(x_n)$ at each step $n = 0, 1, \\ldots$.\n",
    "\n",
    "*Solution*:\n",
    "\n",
    "```\n",
    "n x                       f(x)       f'(x)\n",
    "0 1.5000000000000000e+00 -5.625e+00  7.750e+00\n",
    "1 2.2258064516129030e+00  3.807e+00  1.877e+01\n",
    "2 2.0229637924064581e+00  3.487e-01  1.537e+01\n",
    "3 2.0002760690096917e+00  4.142e-03  1.500e+01\n",
    "4 2.0000000406383567e+00  6.096e-07  1.500e+01\n",
    "5 2.0000000000000009e+00  1.421e-14  1.500e+01\n",
    "```\n",
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exercise above illustrates an important feature of Newton's method. Once Newton's method \"gets going,\" the number of correct digits usually *doubles* at each step.  For example, step 3 got 3 digits right, while step 4 got 7, and step 5 got 15.  That's pretty quick convergence.  \n",
    "\n",
    "Such convergence is called **quadratic convergence** or **second-order convergence**.  In other words, if the error is $\\epsilon_n$ at step $n$, the error at step $n+1$ will be $\\epsilon_{n+1} = a \\epsilon_n^2$ for some constant $a$.  More generally, the **order of convergence** of a method is the value of $p$ for which $\\epsilon_{n+1} = c\\epsilon_{n}^p$ for some constant $c$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr />\n",
    "\n",
    "**Exercise**:  Use the last three values of $x$ in the last exercise to estimate the order of convergence of Newton's method.\n",
    "\n",
    "*Solution*:  Using logarithms, we have $\\log \\epsilon_4 = p \\log \\epsilon_3 + \\log c$ and $\\log \\epsilon_5 = p \\log \\epsilon_4 + \\log c$.  Those are two equations for the two unknowns $p$ and $\\log c$.  We can set up and solve the  $2\\times 2$ system via:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T18:55:46.845285Z",
     "iopub.status.busy": "2026-07-05T18:55:46.845188Z",
     "iopub.status.idle": "2026-07-05T18:55:46.848022Z",
     "shell.execute_reply": "2026-07-05T18:55:46.847577Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.9990273906811293\n"
     ]
    }
   ],
   "source": [
    "e_3 = 2.0002760690096917e+00-2\n",
    "e_4 = 2.0000000406383567e+00-2\n",
    "e_5 = 2.0000000000000009e+00-2\n",
    "import numpy as np\n",
    "A = np.array([[np.log(e_3), 1], [np.log(e_4), 1]])\n",
    "b = np.array([np.log(e_4), np.log(e_5)])\n",
    "p = np.linalg.solve(A, b)[0]\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "That's pretty darn close to 2, and Newton's method is therefore confirmed to exhibit second-order convergence (for this problem).\n",
    "\n",
    "<hr />\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Knowledge of a method's order of convergence is useful when verifying its implementation.  If we saw less than second-order convergence for the exercise above, then the implementation might be suspect.\n",
    "\n",
    "### Secant Method\n",
    "\n",
    "If we don't know the derivative of $f(x)$, then Newton's method can't be used.  Of course, the bisection method is one option, but we can adapt Newton by leveraging the finite-difference methods from [Lesson 19](L19_NumericalDifferentiation.ipynb).  The result is the **secant method**.  The algorithm is nearly identical to Newton's method, but $f'(x)$ is approximated when needed.  The classical secant method begins with two values of $x$ and the corresponding function values.  The finite\n",
    "\n",
    "```\n",
    "'''Secant method for finding the root f given x_0'''\n",
    "Input: f, x_0, tau\n",
    "Set x = x_0\n",
    "Set x_0 = x - tau\n",
    "While |f(x)| > tau do\n",
    "   # Approximate the derivative\n",
    "   Set fp = (f(x)-f(x_0))/(x-x_0)\n",
    "   # Compute the Newton \"step\"\n",
    "   Set Delta = -f(x)/fp\n",
    "   # Store old x, and compute next value\n",
    "   Set x_0 = x\n",
    "   Set x = x + Delta\n",
    "Output: x\n",
    "```\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**: Apply the secant method *manually* to solve $f(x) = x^3 + 2x^2 - 5x - 6 = 0$ to within $\\tau = 10^{-5}$. Use $x_0 = 3/2$. Write out $n$, $x_n$, $f(x_n)$, $f'(x_n)$ at each step $n = 0, 1, \\ldots$.  How does it compare to Newton's method?\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**: Implement the secant method as a Python function.\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**: Implement a function `newton(f, x_0, fp=None, tau=1e-8, delta=1e-8)` that applies Newton's method when `fp` is provided and the secant method if `fp` is not provided.  Here, `delta` is to be the step used in the finite-difference approximation to $f'(x)$.\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**:  Use the secant method *without* defining $f'(x)$ to solve the following equations:\n",
    "\n",
    "  1. $\\tan(x)/x^2$ for $x_0 = 1$.\n",
    "  2. $s^2 = 3$ for $s_0=1.5$. (What is this doing?)\n",
    "  3. $2\\cosh(x/4)-x = 0$ with $x_0 = 4$.\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**:  Another method that can be used when the derivative is not available is [Steffensen's method](https://en.wikipedia.org/wiki/Steffensen%27s_method).  Steffensen's method produces a sequence of approximate roots according to $x_{n+1} = \\frac{f(x_n + f(x_n)) - f(x_n)}{f(x_n)}$.  Try this method on nonlinear equations from the last exercise.  Does it converge quadratically like Newton's method?\n",
    "\n",
    "<hr />\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derivative Information from Lesson 19\n",
    "\n",
    "Newton's method needs a derivative. There are now three practical cases:\n",
    "\n",
    "1. an analytic derivative is available and cheap,\n",
    "2. a finite-difference derivative is good enough,\n",
    "3. a complex-step derivative is available because the callable function supports complex inputs.\n",
    "\n",
    "The iteration is the same in each case:\n",
    "\n",
    "$$\n",
    "  x_{n+1}=x_n-\\frac{f(x_n)}{f'(x_n)}.\n",
    "$$\n",
    "\n",
    "Only the way we estimate $f'(x_n)$ changes. This is a useful design pattern: keep the solver logic separate from the derivative approximation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T18:55:46.849098Z",
     "iopub.status.busy": "2026-07-05T18:55:46.849008Z",
     "iopub.status.idle": "2026-07-05T18:55:46.852395Z",
     "shell.execute_reply": "2026-07-05T18:55:46.851920Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact              root = 0.739085133215161, iterations = 5\n",
      "central difference root = 0.739085133215161, iterations = 5\n",
      "complex step       root = 0.739085133215161, iterations = 5\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def newton_scalar(f, dfdx, x0, tol=1e-12, maxit=20):\n",
    "    x = float(x0)\n",
    "    history = []\n",
    "    for iteration in range(maxit):\n",
    "        fx = f(x)\n",
    "        history.append((iteration, x, fx))\n",
    "        if abs(fx) < tol:\n",
    "            break\n",
    "        x = x - fx / dfdx(f, x)\n",
    "    return x, history\n",
    "\n",
    "\n",
    "def central_derivative(f, x, Delta=1e-5):\n",
    "    return (f(x + Delta) - f(x - Delta)) / (2 * Delta)\n",
    "\n",
    "\n",
    "def complex_step_derivative(f, x, Delta=1e-30):\n",
    "    return np.imag(f(x + 1j * Delta)) / Delta\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    return np.cos(x) - x\n",
    "\n",
    "\n",
    "def exact_derivative(f_unused, x):\n",
    "    return -np.sin(x) - 1\n",
    "\n",
    "for name, derivative in [\n",
    "    ('exact', exact_derivative),\n",
    "    ('central difference', central_derivative),\n",
    "    ('complex step', complex_step_derivative),\n",
    "]:\n",
    "    root, history = newton_scalar(f, derivative, x0=0.5)\n",
    "    print(f'{name:18s} root = {root:.15f}, iterations = {len(history)}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving Systems of Nonlinear Equations Numerically\n",
    "\n",
    "The solution of nonlinear systems is considerably more challenging than single-variable equations.  However, the Newton (and secant) methods described so far are up to the challenge.  A system of nonlinear equations can be written generically as $\\mathbf{f}(\\mathbf{x}) = \\mathbf{0}$, where the bolded names indicate vector quantities.  Whereas Newton's method for single-variables equations requires the derivative $f'(x)$, Newton's method for systems of equations requires the **Jacobian matrix**.  For a system of $n$ unknowns, the Jacobian $\\mathbf{J}$ is a square matrix defined as\n",
    "\n",
    "$$\n",
    "  \\mathbf{J}(\\mathbf{x}) = \\begin{bmatrix}\n",
    "    \\frac{\\partial f_0}{\\partial x_0} & \\frac{\\partial f_0}{\\partial x_1} & \\ldots & \\frac{\\partial f_{0}}{\\partial x_{n-1}} \\\\\n",
    "    \\frac{\\partial f_1}{\\partial x_0} & \\frac{\\partial f_1}{\\partial x_1} & \\ldots & \\frac{\\partial f_{1}}{\\partial x_{n-1}} \\\\\n",
    "         &     &  \\ddots &  \\\\\n",
    "    \\frac{\\partial f_{n-1}}{\\partial x_{0}} & \\frac{\\partial f_{n-1}}{\\partial x_{1}} & \\ldots & \\frac{\\partial f_{n-1}}{\\partial x_{n-1}} \n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where $x_i$ indicates the $i$th unknown (or element $i$ of $\\mathbf{x}$).\n",
    "\n",
    "Then, Newton's method leads to the sequence\n",
    "\n",
    "$$ \n",
    "  \\mathbf{x}_{n+1} = \\mathbf{x}_{n} - \\mathbf{J}^{-1}(\\mathbf{x}_n) \\mathbf{f}(\\mathbf{x}_n) \\, .\n",
    "$$\n",
    "\n",
    "Here, $\\mathbf{f}(\\mathbf{x}_n)$ and $\\mathbf{J}(\\mathbf{x}_n)$ indicate the function and its Jacobian are evaluated at $\\mathbf{x}_n$. \n",
    "\n",
    "Note that if we let $\\mathbf{\\Delta}_n = \\mathbf{J}^{-1}(\\mathbf{x}_n)\\mathbf{f}(\\mathbf{x}_n)$, then it becomes clear that $\\mathbf{J}(\\mathbf{x}_n)\\mathbf{\\Delta}_n = \\mathbf{f}(\\mathbf{x}_n)$.  In other words, each Newton step for a system of *nonlinear* equations requires the solution of a *linear* system of equations where the Newton step $\\mathbf{\\Delta}_n$ is the unknown.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr />\n",
    "\n",
    "**Exercise**:  Consider the nonlinear system of equations $x-y = 4$ and $x^2 + y = -3$.  Derive the Jacobian matrix, and apply one step of Newton's method using $x_0=1/2$ and $y_0=-3$. \n",
    "\n",
    "*Partial Solution*:  Here, we have \n",
    "\n",
    "$$\n",
    "\\mathbf{f}(\\mathbf{x}) = \\left[ \\begin{array}{c}\n",
    "  x-y - 4\\\\\n",
    "  x^2 + y + 3\n",
    "\\end{array}\\right ] = \n",
    "\\left[ \\begin{array}{c}\n",
    "  0\\\\\n",
    "  0\n",
    "\\end{array}\\right ] \\, ,\n",
    "$$\n",
    "\n",
    "where $\\mathbf{x} = [x, y]^T$.\n",
    "\n",
    "The Jacobian is, therefore, defined as\n",
    "\n",
    "$$\n",
    "\\mathbf{J}(\\mathbf{x}) =\n",
    "\\left[ \\begin{array}{cc}\n",
    "  1 & -1 \\\\\n",
    " 2x & 1 \\\\\n",
    "\\end{array}\\right ] \\, .\n",
    "$$\n",
    "\n",
    "With $\\mathbf{x}_0 = [1/2, -3]^T$, Newton gives $\\mathbf{x}_1 = \\mathbf{x}_0 - \\mathbf{J}^{-1}(\\mathbf{x})\\mathbf{f}(\\mathbf{x}_0)$.  \n",
    "\n",
    "\n",
    "<hr />\n",
    "\n",
    "**Exercise**:  What's the Jacobian matrix for the *linear* system $ax+by = 1$ and $cx+dy=2$?  How many Newton iterations are required to solve the system for *any* initial values $x_0$ and $y_0$?  \n",
    "\n",
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jacobian Rows Are Gradients\n",
    "\n",
    "For a nonlinear system, each residual component is a scalar field. If\n",
    "\n",
    "$$\n",
    "\\mathbf{f}(x,y)=\\begin{bmatrix}f_1(x,y)\\\\f_2(x,y)\\end{bmatrix},\n",
    "$$\n",
    "\n",
    "then the Jacobian is\n",
    "\n",
    "$$\n",
    "\\mathbf{J}(x,y)=\n",
    "\\begin{bmatrix}\n",
    "\\nabla f_1(x,y) \\\\\n",
    "\\nabla f_2(x,y)\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "\\partial f_1/\\partial x & \\partial f_1/\\partial y \\\\\n",
    "\\partial f_2/\\partial x & \\partial f_2/\\partial y\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "So the gradient from Lesson 19 is not a side topic. It is the row-level building block of system Newton methods.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `scipy.optimize.fsolve`\n",
    "\n",
    "Although implementation of a Newton (or secant) method for systems would be possible using NumPy, it is more practical to use the tools included in SciPy.  In particular, the `scipy.optimize.fsolve` function implements a method that combines Newton's method with some other techniques applied when Newton's method cannot reliably converge (e.g., what happens when $f'(x)$ vanishes near the root of interest?).  \n",
    "\n",
    "The signature (from `help(fsolve)`) is\n",
    "```python\n",
    "fsolve(func, x0, args=(), fprime=None, full_output=0, col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, epsfcn=None, factor=100, diag=None)\n",
    "```\n",
    "The arguments of interest here are `func`, `x0`, `args`, and `fprime`; refer to the documentation for more on the the other arguments.  From `help(fsolve)`:\n",
    "```\n",
    "func : callable ``f(x, *args)``\n",
    "    A function that takes at least one (possibly vector) argument.\n",
    "x0 : ndarray\n",
    "    The starting estimate for the roots of ``func(x) = 0``.\n",
    "args : tuple, optional\n",
    "    Any extra arguments to `func`.\n",
    "fprime : callable(x), optional\n",
    "    A function to compute the Jacobian of `func` with derivatives\n",
    "```\n",
    "Hence, `func` represents our nonlinear system function $\\mathbf{f}(\\mathbf{x})$, and `args` represents any values needed to define what $\\mathbf{f}(\\mathbf{x})$ does.  The initial guess is `x0`.  If available, `fprime` should compute $\\mathbf{J}(\\mathbf{x})$; if not provided, `fsolve` approximates the Jacobian using a forward-difference approximation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr />\n",
    "\n",
    "**Exercise**: Solve the nonlinear system of equations $x-y = 4$ and $x^2 + y = -3$ using `fsolve`  $x_0=1/2$ and $y_0=-3$.\n",
    "\n",
    "*Solution*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T18:55:46.853387Z",
     "iopub.status.busy": "2026-07-05T18:55:46.853298Z",
     "iopub.status.idle": "2026-07-05T18:55:47.013148Z",
     "shell.execute_reply": "2026-07-05T18:55:47.012680Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6180339887498932 -3.3819660112501064\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import fsolve\n",
    "\n",
    "def f(z):\n",
    "    x, y = z # unpack z\n",
    "    return [x-y-4, x**2+y+3]\n",
    "\n",
    "z0 = [1/2, -3]\n",
    "x, y = fsolve(f, z0)\n",
    "print(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Further Reading\n",
    "\n",
    "Check out the [SciPy documentation](https://docs.scipy.org/doc/scipy/reference/optimize.html) on optimization and root finding."
   ]
  }
 ],
 "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.13.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}