Spectral-eigenvalues.ipynb

{
 "metadata": {
  "name": "Spectral-eigenvalues"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Spectral solution of eigenvalue problems "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Developed by Randy LeVeque for a course on Approximation Theory and Spectral Methods at the University of Washington.\n",
      "\n",
      "See <http://faculty.washington.edu/rjl/classes/am590a2013/codes.html> for more IPython Notebook examples."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This Notebook presents some examples from Chapter 9 of \"Spectral Methods in Matlab\" by L. N. Trefethen. See <http://people.maths.ox.ac.uk/trefethen/spectral.html> for other m-files.\n",
      "\n",
      "p22.m is also rewritten using Chebfun to avoid solving ill-conditioned Vandermonde matrices."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Using Matlab and Chebfun from IPython..."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First we need to set things up to use Matlab and set the path to find chebfun.  \n",
      "\n",
      "You must have pymatbridge installed and chebfun available for this to work."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pymatbridge\n",
      "ip = get_ipython()\n",
      "pymatbridge.load_ipython_extension(ip)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Chebfun is required -- to run this notebook you will need to download and add the correct directory to the path here, along with the directory of programs from \"Spectral Methods in Matlab\""
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "path(path,'/Users/rjl/chebfun/chebfun_v4.2.2889/chebfun')\n",
      "path(path,'/Users/rjl/SMM/all_M_files')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Program 21 -- eigenvalues of Mathieu operator $-u_{xx} + 2qcos(2x)u$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "% p21.m - eigenvalues of Mathieu operator -u_xx + 2qcos(2x)u\n",
      "%         (compare p8.m and p. 724 of Abramowitz & Stegun)\n",
      "\n",
      "  N = 42; h = 2*pi/N; x = h*(1:N);\n",
      "  D2 = toeplitz([-pi^2/(3*h^2)-1/6 ...\n",
      "                 -.5*(-1).^(1:N-1)./sin(h*(1:N-1)/2).^2]);\n",
      "  qq = 0:.2:15; data = [];\n",
      "  for q = qq;\n",
      "    e = sort(eig(-D2 + 2*q*diag(cos(2*x))))';\n",
      "    data = [data; e(1:11)];\n",
      "  end\n",
      "  clf, subplot(1,2,1)\n",
      "  set(gca,'colororder',[0 0 1],'linestyleorder','-|--'), hold on\n",
      "  plot(qq,data), xlabel q, ylabel \\lambda\n",
      "  axis([0 15 -24 32]), set(gca,'ytick',-24:4:32)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAR90lEQVR4nO3d2XKj\nSBBAUZjo//9l5oE2jdmEWCuzznmYcC/2yFIH11mUoO26rgGAaP57+wEAwBECBkBIAgZASAIGQEgC\nBkBIAgZASAIGQEgCBkBIAgZASAIGQEgCBkBIAgZASAIGQEgCBkBIAgZASAIGQEgCBkBIAgZASAIG\nQEgCBkBIAgZASAIGQEgCBkBIAgZASAIGQEgCBkBIAgZASAIGQEh/3n4Ax7Vt23/Qdd38lwDkFnsC\n67qu67pxusa/BCCxwAGbTFoGL4CqBF5CbEbLhuPfGZds8hdEjm+1beNfDZQpdsDGZ7/6/84TJVoc\n0//w458PFCtqwCaTVk+ruIR0QQhRAzbZuzEewhol45Bhvdk/HwghasCa35VSLM4wckFEgQMG50kX\nxCVgVEq6IDoBozrSBTkIGBWRLshEwMjP9kJIScDIzMgFiQkYOUkXpCdgZCNdUAkBIwknuqA2AkZ4\nRi6ok4ARmHRBzQSMeKwWAo2AEYuRCxgIGAEYuYA5AaNoRi5gjYBRIiMX8JGAURYjF7BT4IC1Pz+l\nD7dj7n/H3ZkjMnIB3wocsOanVW3bdl3X/3f45dsPjV10CzgscMBUKi7dAs4LHLDeeN5aXEIcVhrn\nf8TzAp3iUlkoXOCAjXM1ydg4VKJVgjJjMPrZpmlmj62ohwrMBQ5YI07Fe71baw8g0CAIrIkasH78\nGm9EnG9K5BXjsaa0bg38G4EEogZsXindetdj3Zqs+y3+H/1bgBpEDRgluC9aG19ZnICegPGdWyet\n18+ZAYEIGLuc75alP+BaAsaqY9Fa2+AnTsC1BIx/tt8XdccnAhwmYLXb05756p+lP+B1Alad7WK1\n7effASiBgCW3PTy17ecJTL2AMglYKos7/cYs/QFpCFhUH3elu9wfkJuAle7jUNWsVEq6gNwE7H2H\nEwVQMwG70Z4yNeJUKte1gsIJ2C47UzTn2BeObkEU6QPWHW7P9As5nKWmWxBO+oC17hPGBt2CuNIH\nDJZ5mwFEFzhg7c8Pz+MZq22NXGwxckEagQPW/KRriFZ71fkuMjJyQTKBAzaZtPqMaRgTRi7IKnDA\neh/XDMdJs7pYFSMX5Bb4jFFfpsXFw+GbckqsQkYuqETsCWyIk2LRGLmgMlED1o9cixsRqY2RC+oU\nNWBrxVKyqhi5oGZRA0blpAsQMCKxWggMBIwYjFzAhIBROukCFgkY5ZIuYIOAURwnuoA9BIyCGLmA\n/QSMIkgX8C0B42XSBRwjYLxGuoAzBIwXSBdwnoDxHNsLgQsJGE8wcgGXEzDuJV3ATQSMu0gXcCsB\n43rSBTxAwLiSdAGPETCuIV3AwwIHrP3ZlN113fyXPMPOeOAtgQPW/E7X0K22bTXsAUYu4F2BA7az\nUsNktv9T2CZdQAkCB6zXz1vjSk2I1oWkCyhH4IBNVg65lXQBpQkcsEa97mePBlCsqAHrx6/xzkO7\nEK9l5AIKFzVg80rp1lWkCwghasC4nNVCIBYBw8gFhCRg9TJyAaEJWI2MXEACAlYRIxeQiYBVwcgF\n5CNgmRm5gMQELCHdAmogYHnoFlAVAQtPt4A6CVhUugVUTsAiGd/17K1uLW5oXLwdW1F/R+Yhnzb3\nNXD7212+/ShOuTtaH4/1k78Q9Ols2yMhBEpmAivRhc34+KU+fvEcx/QD33jj/XNQNgErxbFJa+Oz\nHHzP8+xByQTsNXvGrLWdGnt2cDj4ArkJ2EPWTjVtn5tZi5A4LVo8rXWeZxvKFD5gwzaN9ufo9e6u\njbVj6JCryd88dm6mZtuV8uxBPQIHrB0dyfqPh5Ld2rCvfsyfPBCH16/YKAhsCBywyeDV/M7YmqtW\nmRxGL+fNW8BXAgdsYujW7wmsc1gsVo53mAFvyROwFeHfyJxMCRcTAXJIHzDeJ1rAHZIErOu6QnYh\n0hMt4G7hAzbkSrdK4Br5wGPCB4zXGbaAVwgYB+kW8C4B4zsWCYFCCBi76BZQGgHjA7dlAcokYCwz\ncgGFEzB+0S0gCgHjL0uFQCwCVjsjFxCUgNXLyAWEJmDVMXIBOQhYRYxcQCYCVgXpAvIRsMysFgKJ\nCVhORi4gPQHLRrqASghYHtIFVCV8wNq27e/F3P6c8Knt1sxOdAF1ChywdnxHxaZpftI1JC09IxdQ\ns8ABq3nwki6AwAFbNB+/xoNagshJF0AvT8D6UM0TlSBajRNdADN5AtZkadWEkQtgUZKA9eNXsvNh\n0gWwIXzA+lblKNZAugA+Ch+wTJzoAthPwIpg5AL4loC9ycgFcJiAvcPIBXCSgD3KyAVwFQF7iJEL\n4FoCdi8jF8BNBOwWugVwNwG7km4BPEbALqBbAM8TsIPGd9PULYDnCdhek/s/ixbAuwRs2SRXze9i\ntW3TttOGzT+lmXXuqr/T2JcPVG96/+Jkhhs0L1Zhw8eoFPK0baxkWuQEcksfsH8fp/5GP7DNBMgn\n/RJi8kLvtPEcFDtcAmxLHzA++Jgr0xtQpmwBa9u2SXeD5ncNz6X1WKAoqQI22rJh5fB6e9YhPevA\nY1IFjEXf7sBcsxEnUxrwvGwBmy8htqNjau6xbC1UV33TH0PYdaIFPCdPwMbLhuOPs0br+d2D21+/\nf3P3xme5lAlwrTwBS6/wAKw9nrXQupIIcFKegHVdN6wWppm6EpxSWrsI1vDxZEQL+m0Cz8sTsCZj\nt7J8Q79sX6YrQbOBZ6QKWGjVHrgnV0kem18xGWAgYC+rtluL5sPZ/PlRNaAnYO/IvUh4lbXFRs8e\n0AjY82xVOGy+I9+TCTUTsIcYGi60Npl5bqEqAnYvx9YHeK801EnA7mJ163mTDY1+eoDcBOxiDpqF\nGI9lXhRIScAuY+Qq02LJGq8UxCdgF5CuECwwQjICdpwjYFzzscyLCOEI2BEOeWnMdzB6WSEKAfuO\ndGXlptIQjoDt4sfzeswvLuxFhzIJ2AcOYTXzukPJBGyZkQugcAI2ZeQCCCFhwNq2PXBrZiMXQCxl\nBexYeyZf4cu///cD3QKIpZSA9eE5X6+u67Yb5oLlADm8GbBxaU6ma103z9n8f+Uu9QDhnF2yO/4/\n/r1aeH4Cmwxe/ZfauSa5J3IAFOW1gDVLdTl/DmzyRS75go3CAZTnzSXEeVperOk2q44ApSllE8eF\nnqngYtLW/giAyyUM2FtcDRbgSQJ2PdECeICAPcoyI8BVBOxRlhkBriJgp+y5cNVin0QL4CQB22ux\nVXs6tBa5xc91LXyAnQRs1VVXTVz7xMWv3//XqTKAjwRs6rF4TL7+Ys8MZABrBKxpythSsd0zACZq\nD1ixi3Xj/YqTB+kqVgBNtQErtltz40c4POxAjx/gJtUFLPRZpfnbyEJ/OwBnVBSwTMf6ScksKgIV\nqiJgmdI1MZ/Juk7PgCrkD1glR/P5TNYkbTZAL3/A+omkqeZovrF9ESCT/95+AE/our8Zq+rNVf13\n3avqGwcqkWoCa3+O04s3Za7zTVSThlX1vQO55QlY27ZDt8YfT1R7BJ9cZbGp+KkAcsgTMPaYnyGT\nMSCo9AHrtmeOag/iw0BW7TMARJc+YL/WEudXfK/8ou8yBsSVPmC/LF7xXcZs9AAiyhOwruu2dyEu\nfUrTzK5hUTO30wQCyROwZne3Zp/194PK6zWwrgiEkCpgJ1W+kDghY0DhBGxqkrHa3vg84fQYUCwB\nWzbOmGN3YzwFylPFtRAPcznBiTqvKgmUScA+m+yzR8aAEgjYXho2IWPAu5wD+4LzQHOeE+AtAva1\nOm/Lsk3GgOcJ2EG13eh5DxkDniRgxzleL/K0AM+wieMsmzsW2eIB3M0EdgFnxdaYxoD7CNhlnBVb\nI2PAHQTsSkaxDTIGXEvArje+r5gj9YSMAVexieMWrj61zRYP4DwBu5GGbZMx4IyEAWtLOhwOx2jW\nyDxwTLaAFVWvgYZtM4oBB6QKWNu2XakbAxygP5Ix4Cv5dyGOZ7J382aT/R62KQI7lTuybJssFXZd\nN/+dptSZzNF5J08UsKHE4/tJ42iVGbCeUWwnGQMWpToHFoudHTvZyQksShiwYkeuOXsW9vNcARP5\nN3EUzs6O/ezvAMYSTmARWSLbz257oGcCK4W7sXzFNAaYwApit8K3PGNQMwErjiPyt6woQp0sIZbI\ncuK3rChChQSsUHYnHiBjUBVLiEWznHiAbYpQCQErnWPxMfZ3QHoCFoBj8WHyD4kJWBgadowVRchK\nwCJxID7MFAv5CFgwDsRn+AkAMhGwkDTsMCuKkIaARaVhZ8gYJCBggTkEnzR+4zMQjitxxOaCHSe5\neAfEZQLLwHLiSVYUIaJsAWvbtq3yIOT4e54dnhBLqoC1bdt1Xdd11TbM8fc8PwpAFPnPgY1j1lVw\nlqM//lbwjd7IiTEIIVvA+lyNQ1VDtCY07BI2yEDhogZsskjYLxsOrRp/XCe3xLyKJxCKFTVglfdp\nDwMEkFvUgM2N927I28ByIpBVnoA1urXCciKQUqpt9KxxzSQgHwGrhXeJAckIWF00DEhDwKrjShNA\nDgJWI8uJQAICVi8NA0ITsKppGBCXgNXOKTEgKAHDKTEgJAHjLw0DYhEw/tEwIBAB4xcNA6IQMKZs\n6wBCEDAW2NYBlE/AWGUUA0omYGxxHxagWALGBxoGlCl8wNqfI2v7493Hk5JTYkCB/rz9AI6bt6rr\nuv73+w+4Vt8wTy1QiMATWNd141CJ1gPMYUA5Ak9gi+bj13hQE7nzhoZ5LoF3xQjYZLVwsUP935n/\nkWhdbtjW4akFXhQjYDsjpFVPckoMeFeMgH3Uj1/DoKZkz9Aw4EXhA9a3SrHe4pQY8JbAuxAphHeJ\nAa8QMK6hYcDDBIzLaBjwJAHjShoGPEbAuJibsADPEDCu5wL2wAMEjFvYmgjcTcC4kYYB9xEw7qVh\nwE0EjNvZ1gHcQcB4glNiwOUEjOdoGHAhAeNRGgZcRcB4mlNiwCUEjBd4pzNwnoDxDts6gJMEjDdp\nGHBY+IC1P8e/9se7j4dvaRhwzJ+3H8Bx41b1H3dd13/cub99KEPDvG7AfoEDNuRq+J1xxohl2Nbh\n1QN2ChywiaFbkwlsXDhtK1w/inmVgD1iBGxyZuurDolWLBoG7BQjYCJUFQ0D9ogRsI+6rhumNLVL\nwLYO4KPwARtypVvJ2NYBbAv/PjBy8y4xYI2AUToNAxYJGAFoGDAnYMTgJizAhIARhgvYA2MCRjAa\nBvQEjHg0DGgEjKA0DBAworKtAyonYARmWwfUTMAIT8OgTgJGBhoGFRIwktAwqI2AkYeGQVUEjFRs\nTYR6CBjZ2JoIlRAwctIwSC/wHZnbn+NTfy/myS+hb5h/DpBV4IA1v9M1dKttWw2jN8xh/kVAPoGX\nEFWKPfp/JpYTIZ/YE1jzM2+168en8R9pXp2Ghnn9IZMYAZv0aXHlcI1o0XNKDJKJEbC1CIkTX9Ew\nyCRGwOb68Wu889AuRPbQMEgjasDmldItdtIwyCHwLkQ4zBWnIAEBo1KuOAXRCRhV0zCIS8ConYZB\nUAIGGgYhCRg0jW0dEJCAwV+2dUAsAga/aBhEIWAwpWEQgoDBAg2D8gkYLNMwKJyAwSoNg5IJGGxx\nzV8oloABEJKAARCSgAEQkoABEFLggLU/Jr/51uO5SfTvyOMHbhI4YE3TdF3Xdd1wiHGsAahH4IB1\nvzc4t23b2fIMUI0/bz+Asz52K8FYFv1biP74gTLFmFomR8D+Mfe/Of548hcASCxGwBYtzl4WEgEq\nEXUJsR+5hsFLtABqY14BIKTAuxABqFnUJcQ9oi8wjnemRPwWhvORQV+I8fnU6K8FpJQ2YJOjT9CD\nTtCHPT7cR3whFvf9h3jkUBVLiEWbXysrhP4KKW8/iuMWH3/Q1wISE7CiTa6VxYu8FlAaAStX6CEm\nGa8FFEjACuUn/XJ4LaBMaTdxjFd7Iv74HP3xDxJ8Iwm+BUgpxq4wAJiwhAhASAIGQEgCBkBIAgZA\nSAIGQEgCBkBIad8HxkeusA6EJmCVWrtXCEAUlhABCEnAAAjJEmKl3BkEiE7A6uUcGBCaJUQAQnI1\negBCMoEBEJKAARCSgAEQkoABEJKAARCSgAEQkoABEJKAARCSgAEQkoABEJKAARCSgAEQkoABEJKA\nARCSgAEQkoABEJKAARCSgAEQkoABEJKAARDS/2uDXzbUsmXQAAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Program 22 - 5th eigenvector of Airy equation $u_{xx} = \\lambda x u$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "\n",
      "% p22.m - 5th eigenvector of Airy equation u_xx = lambda*x*u\n",
      "\n",
      "  for N = 12:12:48\n",
      "    [D,x] = cheb(N); D2 = D^2; D2 = D2(2:N,2:N);\n",
      "    [V,Lam] = eig(D2,diag(x(2:N)));      % generalized ev problem\n",
      "    Lam = diag(Lam); ii = find(Lam>0);\n",
      "    V = V(:,ii); Lam = Lam(ii);\n",
      "    [foo,ii] = sort(Lam); ii = ii(5); lambda = Lam(ii);\n",
      "    v = [0;V(:,ii);0]; v = v/v(N/2+1)*airy(0);\n",
      "    xx = -1:.01:1; vv = polyval(polyfit(x,v,N),xx);\n",
      "    subplot(2,2,N/12), plot(xx,vv), grid on\n",
      "    title(sprintf('N = %d     eig = %15.10f',N,lambda))\n",
      "  end\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "text": [
        "[\bWarning: Polynomial is badly conditioned. Add points with distinct X values,\n",
        "reduce the degree of the polynomial, or try centering and scaling as described\n",
        "in HELP POLYFIT.]\b \n",
        "[\b> In polyfit at 76\n",
        "  In web_eval at 44\n",
        "  In webserver at 241]\b \n",
        "[\bWarning: Polynomial is badly conditioned. Add points with distinct X values,\n",
        "reduce the degree of the polynomial, or try centering and scaling as described\n",
        "in HELP POLYFIT.]\b \n",
        "[\b> In polyfit at 76\n",
        "  In web_eval at 44\n",
        "  In webserver at 241]\b \n"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO3d4bKj\nqhKGYTy17/+Wc34wcRhAgojwNbxP7dq1ZpmYXtLagEaPz+fjAACw5n+zAwAAoAUFDABgEgUMAGAS\nBQwAYBIFDABgEgUMAGASBQwAYBIFDABgEgUMAGASBQwAYBIFDABgEgUMAGASBQwAYBIFDABgEgUM\nAGASBQwAYBIFDABg0n+zA7DhOI7z0dXhzw3v9T88fxD23TDqo0qDvAq7/Mrz5/S956cX1hCt5OoF\nhZXwtPEu2pK/kDPTk99VpNnVB9UkpMtleLr+dB9s2CXL8S+PAlarYZ+Jctc7M/JhqjW/PYqqXMn8\nvpE9fqU/F45055rDTw932vD10Z9WfkG6kvDYsdv+/JK2LZm2QnaPaNClWctZXQi1JiGzoWb3guhD\nz2oUfUrhQ590r62jgNWKDtCVb3EXHbS7sl2waJRTufKf+6f7tyo8l+51V2Us/U02hvJb0F1z8ofS\njK3UMfnT1V698Wo3+ZmQV28s7wU167wV/yYoYE8VOlyFt9ydhHS5LthV37YhpOzrm3tz9Qep7GeV\n/9L0BVfd58337QF+Zlqvub43kr+mihRCyq7qYQ+1MLoqfGjbJ66BAnZDtv/YXIruaqgKg4WzH+Fs\nyZOee7rFymvbeTrlVXeTP3tS5/x/l9n4qzg7vqzw3lunnQp/cmXGNnzoDihgT90d7rRl3q2zCA9H\nYM26dLfTNURThey6OsqZlu12NLRg9+R/MrWQdrDeTsgpH2oFBeyetB9an0lhJ/TWGz8XVxldnVG4\nO+KJVh79M/uCz+eT/v7qlfUfHf0m3WKVH9qwHfBTffI3p3r5Q3slf1vORCfertLv6o31wbTtkj/j\nXxKVfBF0yrAtkn9bjMBs27bnBZD8oOcCADCJW0kBAExaagqRKQUs7Cq9SXtsa50CdvWFR2AB5e/w\nkvbY0zoFzOMWDNgQaY89LVXAau6/AtwlXhgKIzPgCfHMd4sVsCy1NlCb5yGeMruVQG0zEk+ZWkgm\nMp+rEAEAJq0zAqu5pwtgVP0NvYB9aA1au1MblcMcoylkNGzoMJFCTCECAEyigAEATKKAjaZ2bQ/x\nYIBhzXocf/9TiKeeYEj61rmIA8DO/PE/PGtzHE7+JA4eYQQ2mtp5UeLBAGOaNfqQz+dyHCaYZoIh\n6aOAATCPwdaeKGAAbCtUr8IgDAuggI2mdqqWeDDAxGbN1jDBNBMMSR8FbDS1mW7iwQBqzaoWj5MM\nSR8FDIBhNWe/mEhcFQUMAGASBWw0tZlu4sEA05s1GoRNjyclGJI+ChgAq+5ePU+NWAwFbDS1U7XE\nY8XxdbV0cDy3KDRrGIJCPBHBkPRxKynAgPDZFulzLsSrF/ASRmCAbSae2/SGhrtvcDniYjYqYCKJ\nq9ZZJp5VhfON4VaN5iHHLAr/PzGM9J9zw9BcVJipVrN43+3fiRfnktt9AmUi45urKcToQFOYZlxP\n8/0PuXFiDRMptMs5MJ+yRnoVQC2dikUHEeNtNIUI2PX5fM65HV+opCZ5xo9pnnwifdllUMBGkzru\nOOKx4/N1/jNaOiMo54JaUigMNOtPbKIGFDAAgEkLFrBCR0Zh6kDtvCjx7OON5K+cylNrVoVDQURt\nE5mwWgFjGA5kjbmOaUxh4DJCeEsVsOyFWNF3PrK/ZxGL0kXRC0zjiJ8SHIThLgNX+tfzBezIfWPm\n3IGn78nZKjsR8ZSpxVMpCjtM+767QLq27Pr7bsbnf8L3WCFU19UyTS2erHVGYMfFF/7VqOUE8aBZ\nfQFQa1Yfj9QgTG0TmbDOF5mLtzoV6mcBa/NV4b09jt0Zp3VGYACyOOIXSA3CcNeCBUx8JK42vUk8\nu5lyyFZrVrV4nGRI+hYsYGV0uACEOCbYtV0Bm05tgEg8aHNrZlKtWdXicZIh6aOAAejsvTHNS+fz\nGIQZRQEDVpY94nO8xhooYKOpnaolHiuubg5i4qYhauGl8Uwv6mqbyIR1vgdmhdpMN/GYkL2/zOl8\nQpjs1lMLTC0eJxmSPkZggG1TDnxLfrds+iAMd1HAgBXkbkDzd16x+72Prz7ofJcvBt3vs5y+cvDf\ntcMi/Rnpk+6cQxdH7padczuPavM8xFMmEk9hCvFbMOLqFd3GOlnho72g5u3/3kG4w2bsuOcW4pl1\nfBDJtJNaPFmMwADz1I7FwBhcxDGaWqeGeEzwzwk6f3bBw4NcMAtUv/XevuVu8nFazVqIZ/CWCT5X\naxOZsGMBm5WgwBPZKw856r2BQ4QVTCECAEyigI2mdnkP8eCuhtHJ82bt/SDpH/GMv6SezG9AARtN\nbc6HePCSsAaoNWtNPINrmNomMoECBgAwaf0CxslY7Knyq1oo4N4c4tYvYGrUZrqJZ1sjO3ZqzVof\nz7AapraJTKCAAUAJ4zBZFLDR1E7VEg8GeNis3U8E3I1nQA0j8xtQwADccKuWrDR2WelvWcZSd+Jo\nuKEOAFQ6axgHGBGrjcA+n09417jrl03rTKmdqiUeDKDWrM3x+NL1xl+jtolMWGoEZmLgpRYk8WAA\ntWZ9Ek9Ywzr+WWqbyITVRmDu4sl+0Q/pP2UfLseiWYuiF1jD8fBd7w3FUM/AI8vqHcmT/Y7juNqT\n+YIzaqT9ITXZU7/16X13R3j79W+s4VXB7bKmxtGbfua79UZg+ltcrVNPPKb5o0zlqd+J1GLrGM/n\n83c0dv43N6R9rHMO7Gh9sh+ASg2DoU2uPg83S7aGcUB6wzoFzErFUouTeFYVzqiH00HpTHv2Ldl3\nOXe5qPyun2H0ilB2UUUV/4hVegN7ooFZzic4B4aHxM8ERIfO7M+/1nBjR2jbax7ua+yqU4hnvrfa\nOTAAwCYoYKOpnaolHtP8tRuecn+5uVlfGn4JpplgSPrWOQdmhdpRhnise7jF/EUWb291tWZVi8dJ\nhqRv3xHYJhdHAcCq9i1gfT35/gewNjqLeAkFrA//ZUa/o5b3VbWZbuJBpZWuBhRMM8GQ9FHAOggT\nb5k9HADEUcCeSrul5QkTtVO1xANZ7435BNNMMCR9FLBHrnYwJv0B4G0UMACASRSwduX5jatB2IBT\nteFdsX9eHql26lgtnh1sOGEgmGaCIenji8yjvTrTfVw8JdbXsIvZTq2Zd7V40MWYr0vXE0wzwZD0\nUcBeNHinLXzWSw9Bxz6kyg/gMYXYSG1/ronn/KYaoI9ExU8UsHelBeONme5bq4zqnNrMu1o86KKt\nWd/rIwqmmWBI+rYuYGsMR/zY6+6uvsAfDmBzWxewZrfmD6MyKXKq9oxKJJ6TWjw6zsemVP5eilqz\nqsXjJEPSRwGz7cmpuDUGoJvwj/vy0lp19ft6A5KBfEN3FDDDnl9IwjFlAfTcsS0K2AhhndCb6nk9\nouwXq68+VG/72JA+kTmcVwy3ajTfWLPoOP4svPWucNHVzOf1u/52zh4Gn10UbpkuK1xpkf6M9En6\nMeTPHcdR7qE2DGLaxj3dL7vvu8I3vhVw7gJXa/75AgVpYRjzoeE//fTgGUYUkn9xWr3uhl3IgV7p\nUb8eta+pbGhK5t+11BeZz90+2O4qTaB2J4JXVX5j+vx69UvfsM52Ik00wd2sVclyYKx1Clihi7qe\n7rUwuCjx6aruricsY726+eGaby2VFV6j4XPbJ3k0Qbd22gOR9c+BRfPd2d/fWhStvH6F3+c1X07H\nV67w+QmJcFHw/+MsY60r/FOEPp/bYWSfZ303jO/bj/OLcem7gkU/Puv7ApWScF6FeP4z/GW4SNDZ\nyiLnVgTP8QiGpG+dkUp2BPZzKHa31/9wlPA9vD7a5r3PfqUn/xu/Ft1p/HT4mnG3XRoCyI7Gol8a\nHc23hX3V9M9TLtglf69qwGS7YLOqhaQWT9Y6U4htxp+akqpeLhfPrenE7qevgl5I1ZqfBBBc55b5\nJbwuKad2KFSLx0mGpG/3AjaYlUs5ztm88gWc54tfiqHwKX0/Xb9FAKTWKWDpWe71jCx+V/Vj8GAl\n+7nDPh2AsnUKmHu/bnWqH41zyy9Vr3I0af2YdXKCcmWa2gkVtXicZEj61r8KUc15KltE5T7jL94b\nMuRiH57mvasEg4snf3zEmGkGwTQTDEkfBWwCneEXUEbiQRkFzAAOIgCQooDV6nefCL7RWaIWD7pQ\na1a1eJxkSPooYNNU1jCGXxjp88+TE17JvULmk+24hQI2ejB061TtgP1Z7dSxWjzoQq1Z1eJxkiHp\no4DNVK6d9EYxRc3X2AEFFLDJrmoYhw8AKKOAVelYTtJTtdFd2Ad3ftVOHavFo+N7j/z89um73bp/\n5y8KL9tv2zntnWRI+pa6E4cJV7eZOF57rmNDPBOpxSMi+7CFcOmMoG5Qa1a1eJxkSPooYCrIXrQ5\ngidbAlthChFYUzjf+PAZqm8siv4Vzh+KRLjtosJMtZrFbx+ZTrZcvOzHc0NefYDkXMRTNiWe5IzR\n52oKMX1l+hoF2XjOPWv8JUtq28fphaQWTxZTiM7ZeUwXNlF/4JCtWDV07kcDoyhgo6kdZYjHhPRx\nd7Yq1lWos/4CwU0nGJI+ChhgQ3SAK/8T2AEXcfzA1CIAaKKAjaZ2eQ/xYAC1ZlWLx0mGpI8C9sew\n88lqUz3EgwHUmlUtHicZkj4KGADAJApYCSfAAEDWUgXMxHfI1cIjHgyg1qxq8TjJkPQtVcCcc5/P\np/m+cHytEgAMWaqA9T0L+trz1LUmJYkHA6g1q1o8TjIkfUsVMC+6Q0E4qSh790wWCS7Sn44GNmfp\nbjSh6Mhy3lzH/duROe7fbscPvLh8A15DCikwGjZ0mEghqyOwz7/C3z9e87tnwtQ69cSDAdSaVS0e\nJxmSvnXuheib/0yCJ5Xs1W6HWqeGeDCAWrOqxeMkQ9K3TgGj+QFgK+sUMGBtV7MLXWYdAIusngOz\nS22mm3hM8GfU0685nhcuNX/9cQy12NTicZIh6WMEBpiXXn8L7IAR2GhqRxnisS47MnNi34AMv+gy\nMYzzn2ea6Xz10AXmbqjoBcoMXOn/xGHhqwxQNiWFosOHL07hMffWz0ADEynEFCIgR//AAShgCnE0\ntbE58ZjgB2FedN+Z6Pea1JpVLR4nGZI+6aR/Tnyvhj6jKWQ0bOgwkUKMwAAAJlHAAAAmUcBGU5vp\nJh4MoNasavE4yZD0UcAAACZRwEZTOy9KPBhArVnV4nGSIemjgAEATKKAAQBMooCNpnaqlngwgFqz\nqsXjJEPSRwEbTW2mm3gwgFqzqsXjJEPSRwEDAJhEAQMAmEQBG01tppt4MIBas6rF4yRD0kcBAwCY\nRAEbTe1ULfFgALVmVYvHSYakb8ECJj4SVwuPeKwoPOtd/xnwauGpxeMkQ9K32hOZSQIsKXw4U/Sg\npvARl/TisZWlCpjfgaMaJljS1EIiniWpbUbi+UkwJHFLFbAUHVLswB/4wmwn87EDqwUs6qqcA6/z\n/+zAsCtN78Irr6YWgeVZLWDpjspujGWQwEANqwUM2Ep4cje8ZCP9PbAPBisAAJPWH4HNnVEs9I7D\n8xyzIlTovF/FoLB9zkgs9vM0M1+hWUn7SvqZv3IBm35N6s8T7DrH5VmZWo5h+s4zPYXaTA9buVlJ\n+xrTU6jSgnfiOPkzBLOjKNG/gcJc07ePfgpliYc9vVnFKWwf8RQ6rVzA9EUn4RFh+yyJZi1j+9Rb\nZAqx/nszOjGY6OBMxPapYS7zadYyts8tixQwhVa/FYP+2dG52D6VFLZSfQw0axnb5671t5fatVjh\nrVejReMJxiC1fTyjhxWpzJdqVs0Yolu5KqScfuarxwcAQBYXcQAATKKAAQBMooABAEyigAEATKKA\nAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEA\nTKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyigAEATKKAAQBMooABAEyi\ngAEATKKAAQBM+m92AGYcx/H5fNKfa97of/Bvif7ZK6TmNZw/F8LLflDNK/1rsmsOPzr7muwLwk/5\n+YKrvxq3NGd++vowH3qF1LyG8+cwl6LVph9U88roz/yZ+TW7xlUwm6c9BeyG5vy4OoI/zLYuyXq1\nA2SLRPa957uiF4e/9z+E/0yPa4UDZbjmq8KZvhK9tOXqz3xo1jfzs4lXSKT0r4jyM31BIfPTT8/u\nGlEM6bat/KsXwxTiDT8P6Ffvin4TZW294+v8Z/T7hvDK74p2p2hRtJ7yX3S19Ocbf44AoiPFbp3Q\nARoyv+NoIM3t87DenPauNfPL6Xe1qqsAKjM/DSb6zbZpzwisg8JcVvqCJ7Mx7qL3d9UTLEcV/r5L\nF7scc83rb70M01XmWPR6hbQP11kfTPYjyjFXvqUtElDA7sl2RX/ukM27Sqp+JTWHiS4zOWdPNpwJ\nzB5r0u0QDZ7OPT8cXEbHqfBTXGvpxV23Mr8yH+r1Tfv6lxXeG/2NLqlb2SnEcyXhTpGmvSOx61DA\nOih0+vpm4a2hTM24sDm26L0/j001n1U/FcO+LeIqx/q2VN+0fxJV+sboL00jvDVDjrsoYLelaVqe\npjhfHPXC7mZtuqqreCrXn67w6iNCx/c8c/mVNX94WgXTt7Dz67iV+YX33mqy7mmfXWdNbFGJqkz7\n7FuyEw9XL0ABm2kp5D02RNpvixHYCppHdYBdpD3ouQAATOJ7YAAAkyhgAACT7J0DK0x8MyeOhV2l\nN2mPbRkrYIWvQ159Yx9YQPkmFKQ99mSsgP3UcK8awDrSHntaqoClXVRuL4bnxAtDYWQGPCGe+W6x\nApal1gZq8zzEU2a3EqhtRuIpUwvJROZzFSIAwCRjI7DCfcOab7YG6Lu6USRpj51pDVq7UxuVwxyj\nKWQ0bOgwkUJMIQIATKKALes4nIWzsEBnpP0+KGCjDbi2x5cuP/r/+Wlq1xqpxYMuxjTrmfk/e2+C\naSYYkj5jF3Ggkq9eZw2Tn8oGngrznMzfBCOw0d4+L5rutL5DOiueu9TiQRfj094VM18wzQRD0kcB\nW8pVl7Ncw4BVkflro4Dtgj0Zq2KqcFsUsNHeO1XbthurnTpWiwddTGzWbNdNMM0EQ9JHARtt4kx3\ndk9Wm3lXiwddqDWrWjxOMiR9FLBF0HvDnmomHpg/XxUFbB01HTj2ZADLoICNpjbTTTwY4KVmrT/v\nG3XdBNNMMCR9FLDtMAjDnsj89dgrYMdX4QUj47nrjVO1Ty4jVjt1rBaPjnLmb5j2T6jF4yRD0mfs\nVlLRQ9PTJhffjUX4rij7iyHlzCftsSd7I7ACEw+wAfraNu0bOmHMIi5mqQKWFc66hB3VaDZm2KL0\nZY9X6Jy7/a5zT/a/mbU1lBcV5utMkMr88P8Tw0j/KZVyIosMZb6xvttxPZESbfHwsesjIxyveTKQ\nWcQaIil0lfnZtHcyYb+KzH+ViRQydg6sQGTX9ccT+XZ3jjNhS5BKezcw80ldOHNTiJ/P5xzhXnVC\nJ/I7FfPs6C7NfJ209yofoKqAPXQZxgqYc+7zFf4mesHwoJyr7hL2Pe487IcK7slqx2UdUeaLpL0L\nkrAQAs36E5uogb0CBgCPUS9WsGMBO47Xhx22HgUrhe3znjcyv3IOQK1Z1eJxkiHp266AvXTu1+4p\nZcFZRLxhTIqOSacufwuZv4DtCpi3Ru7arZoYjFTBkvYqYCN346saqTb1fhyHVDlX2z6LmdXWas2q\nFo+TDEnfXgVMgdpMN/GgWbZHaOjB31JdN7VNZMK+Baxj7i4wPyO1J6OvBfIztNifgyf2LWDAnuis\nhNgapm1UwMZ33LL7Rq+Z7l5/zhmPyJ7MmYAlqTWrWjxOMiR9GxUwAB0tM5Un0nVDg60L2JTEVTtV\n++9NuebvyWrbBw3SRFJrVrV4nGRI+rYuYMDylhkneS/9OQpdNzSggD212AGCPXkHtDLWYK+AFR4Y\nauJZol3C61g11TaXWjw6rtJ7n7TvKI1nelFX20QmGHugZeGJzN7Vg5hnjZPSh0aqzXTntuHMYaXa\n9hFRznz954+rBaYWj5MMSZ+9EVgBGYANTUn7xWbOvemDMNy1VAHzon5oOLsSDtK/D7f9k7LpovCf\naovSV/b9rHChyJ88fpH+vFwoN+tQyvzwnw2Lrj7ofJffs/o2zXH8WdhrhQ1/1w6LDGW+7pxD1lGc\nSPnuOZ/09YUO48O+5M+3Ry9Iw37jQ++sKh/PrC52l+3TkUg8hcxP095VZP7baR+9xkra9/2UW0Qy\n7aQWT9ZqI7DBW3zJiRSP6RRD1I7FwBjGLuL4fD7n2DbqkB7f+YVoqRqpfqgrxjPlag7ZhpsrzfyH\naT+4cdWaVS3tnd4mMsFYAXP5q+Y+2d+LmHtRH5aRvfJQNu1NY5+1YrUpxCxysRkTiQBkbVHApKhd\n3vMznsE1TG37IFXZIwwz53mz9u2GqqW9I/ObUMBGZ6ranA/xYAC1Zq2JZ/MjgwkUMPzGRCL2ROaL\no4C1m3Jqbdb5PPZkcyq/qgXYRQEboe/JgL7q4xlTw9S2z8JG9oTUmlUt7Z3eJjKBAoYbGIdhQ6S9\nrPULmNo19Gqnau/G8/bOrLZ90MXDZu2+F6ulvSPzm6xfwPAGOqTbulVLVhq7fG9PPDsOBChgzq21\nmw3w+bDFsCNfucl8HRSw0Z6cqn1jOrQ5npdqGKeyl6TWrE/SnszXQQFr1FxL1Ga6n8TzxqSK2vZB\nF2rN+jAeMl+EvZv5GrXq7UHDSZX1/jrLOB6+i8xXsFQB03+WyqqicwNs/pGepL2JfpVyhNFQTDbO\nVa1TwMoPa9ahFlvHeILH7+Z/PzieHZD2bV5K+zDz765ebROZsE4Bg45wNyyfKmCHtaVhMLTPBavR\nliHtB9ihgB3H8c9DbP/+9mLWJVpUeNetRc79eQJh6wrDvl63CM8t88af/F19ftFFbfuIHe+MHmyO\n7+mZnpl/bo27mf8zjLn75vjMzyHzb1tn0JqdS7k1Kr/Vu2zoij6cylc+E7Aw8YmdqynE+rDfTvvm\nd/V6O9qIZ77HZfQt2KMAYLp1CpifmvM0Ow7fC5a0pgmIxzT9tPeam/WlzqJgmgmGpG+pc2DKO/BJ\nLUjisa7Ld3Lf3upqzaoWj5MMSd86IzAAwFYoYDZw1g127XMlPQajgI2mNtNNPKi0Ui9KMM0EQ9JH\nAQMAmEQB+2vARMfn49S+Hqh26lgtHkz03phPMM0EQ9JHAbttpYkUALCLAmYAJRMAUhSw3amdOlaL\nZwdDrhLUalitaJxzkiHpo4DtTm3mXS0edKHWrGrxOMmQ9FHAAPzGPDYEUcAAKGJGDT9RwP7ByYDp\n1OJBF23N+t6YTzDNBEPSRwG7Z/xEClM3AJBl7270Zz8lPedZWKTj+nmsc6htK7V4dFylN2nfQHBb\nCYakz1gB+/nM2YYHMQP6ypn/PO0HPFRlzHNbsJWlphApWtgQaY9tLVXAvKgf6h9We/6c/X39ovSV\nd1eYru3hCh8uOv8/NwzBRdELxKXDr46Zf+bIrXeFi9It+etdf8dqb7R1uGV0Uk5kkaHMl55qizai\nf3p6YSLFvz6qXnf/wPIsR685kPr1MOsyV0MKdfnQ8J/lzE/TPn1N3YdeZhppv6EpmX+X9Dmwhs33\nfIsXZuq5BBFj3E1j/QMN8AbpApbyXdHzZ/+D7ylEMxXs0lhJmvmkPWCsgLncLup/Y2XXVRuYE48V\n2SsPO26rV68S9M2qcyGiYJoJhqRvwYs4AHTEzDlkUcBGU+tkEQ8GUGtWtXicZEj6KGBz1Nx0kX4o\nABRQwDKy1YVyAgBSKGCjqX1DkHgwQHCd5O+5hwEE00wwJH0UsNHUZrqJB6e0uvSaeKhv1jFTHYJp\nJhiSPgoYAMAkClhe1BV9Y3BfnkvhlBsAlFHAavUqJ2oz3cSD0EvnqNSaVS0eJxmSPgoYgLyXpgEK\nNZKJB9xCAbt07mZ9dyq1U7XEgwHUmlUtHicZkj4K2A+vDuuvuqL0QzGRT0uSEPrs3cx3JHZg7InM\nhwn2RmA/Hxgqfi5ULTzisaKc+eLbLXlE5+Sb3QhuLsGQ9BkbgZWfyOwsJEHyUIx4vx08daM2864W\nj4iaZ5ErU2tWtXicZEj67I3ACnigDjZkNO3Tr1oa/CMw2VIFLCucdQk7qtFszMRF/0b7dzfWiXDP\nRYX5OhNsZf5V5LLBL7zIUOZL992ijegfoH41kZK+OH2NgmxIZ90a3w9V20TE4+5kfjbtnZ3N6BOe\ntHd6IanFkyV9DuzW5pPddWuc0ynWAscr6hPYdNo7mZvTwyjpApbyXdHzZ/+DrV33KtRZf4HaplOL\nR0Sa+WukvZuU+YKbTjAkfcYKmMs1c3JdH3mABZXznLTHhta/iAMAsCQK2Ghql/cQDwZQa1a1eJxk\nSPooYKOpTfUQDwZQa1a1eJxkSPooYAAAkyhgAACTKGCjqc10Ew8GUGtWtXicZEj6KGAAAJMoYKOp\nnaolHgyg1qxq8TjJkPRRwAAAJlHAAAAmUcBGUztVSzwYQK1Z1eJxkiHpo4CNpjbTTTwYQK1Z1eJx\nkiHpo4ABAEyydzf69HEqNYsA667Sm7THtoyNwPwzkLxoytj/M7tIilpsxGPCVeaT9m3U4nGSIemz\nNwIrO/fn2YEA45D22JOxEVjZVRf1/OfV70cuKoQ3JcLz8b5zw0h/nhvGyckrzEnoZP6ZZnPDOP8Z\nPs+dzI8WWcl855z0U8mjjeh30TDzwuCziw5Tj12HoCkpVJ/5NT8DDUykkPQUov7mA95A5gM1jE0h\n+q6oF3Y2rxYJUhubE48JaXqT9k+oxeMkQ9InnfTPie/V0Gc0hYyGDR0mUsjYCAwAAI8CBgAwiQI2\nmtpMN/FgALVmVYvHSYakjwIGADCJAjaa2nlR4sEAas2qFo+TDEkfBQwAYBIFDABgEgVsNLVTtcSD\nAdSaVS0eJxmSPgrYaGoz3cSDAdSaVS0eJxmSPgoYAMAkClQ+TysAAAJJSURBVBgAwCQK2GhqM93E\ngwHUmlUtHicZkj4KGADAJOnngWWd/ZTsOc/zGRNDY7pDLTbisaKQ+aT9XWrxOMmQ9Bkbgfk7/Gcf\noB4uVR6Mq8VGPCYUMp+0b6AWj5MMSZ+xAgYAgGdvCrEsnUsR7NeohUQ81mWnENU2I/H8JBiSOOkC\nFjVneY44fH7o+TPTyrCoPvOzaV9+C7AM6QLGTog9kflAjcPcrpJei3V2PMsXKAKmRelN2gP2ChgA\nAI6rEAEARkmfA+siPLM95dP9D1dfPr1aOobC7NNVDArb54zE4kSFZuYrNCtpX0k/81cuYNOvSb26\nQuykc1yelanlGKbvPNNTqM30sJWblbSvMT2FKq08hehvTzA7ipLjOKwkyhTTt49+CmWJhz29WcUp\nbB/xFDqtXMD06d8BaC62z5Jo1jK2T71FphBvfeVZJAYTHZyJ2D41zGU+zVrG9rllkQKm0Oq3YtA/\nOzoX26eSwlaqj4FmLWP73LX+9lK7Fkvq+6eCMUhtH8/oYUUq86WaVTMGv4kUYjvpZ756fAAAZHER\nBwDAJAoYAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAkChgA\nwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAk\nChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADApP8D\nw4QcY8N6p1MAAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note the Matlab warnings from using polyfit (which uses the ill-conditioned Vandermonde matrix). \n",
      "\n",
      "To avoid this, we rewrite Program 22 to use Chebfun in defining the Chebyshev points and D2 matrix, and to interpolate the resulting eigenfunction and create a smooth plot.  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "\n",
      "  for N = 12:12:48\n",
      "    \n",
      "    % Use Chebfun:\n",
      "    x = chebpts(N+1);\n",
      "    D2op = chebop(@(u) diff(u,2));\n",
      "    D2 = D2op(N+1);\n",
      "    \n",
      "    D2 = D2(2:N,2:N);\n",
      "    [V,Lam] = eig(D2,diag(x(2:N)));      % generalized ev problem\n",
      "    Lam = diag(Lam); ii = find(Lam>0);\n",
      "    V = V(:,ii); Lam = Lam(ii);\n",
      "    [foo,ii] = sort(Lam); ii = ii(5); lambda = Lam(ii);\n",
      "    v = [0;V(:,ii);0]; v = v/v(N/2+1)*airy(0);\n",
      "    \n",
      "    % Use Chebfun to interpolate and plot:\n",
      "    d = domain(-1,1);\n",
      "    vpoly = interp1(x,v,d);\n",
      "    subplot(2,2,N/12), plot(vpoly), grid on\n",
      "    \n",
      "    title(sprintf('N = %d     eig = %15.10f',N,lambda))\n",
      "  end\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO3d7ZKs\nqBKFYerE3P8t1/nBbpvmS1CElfg+MbGnu63SLElNREs/3+/XAQBgzf9WBwAAwBUUMACASRQwAIBJ\nFDAAgEkUMACASRQwAIBJFDAAgEkUMACASRQwAIBJFDAAgEkUMACASRQwAIBJFDAAgEkUMACASRQw\nAIBJFDAAgEkUMACASf+tDsCGz+dzPLo6/PnCe/0P9x+E3RtGe1RpkKWw6688fk7feyy9ModoJqUX\nVGbC08aHuJb8lZxZnvyuIc1KC2pJSJfL8HT+6TZ4YZOsx789ClirC9tMlLvekZE3U+3y26Oo6pXM\nbxvZ/Vf6c2VPd8w5XHq40Yavjz5a/QXpTMJ9x9u254dcW5NpK2S3iAuGNGs9qyuhtiRkNtTsVhAt\n9KhG0VIqC73TvbaOAtYq2kE3vsUVOmi9sl2w6Cincean26f7WxXuS7e6UhlL/5KNof4WDHc5+UNp\nxjYamPzpbEtvLG0mpwlZemN9K2iZZ1f8L0EBu6vS4aq8pXcQ0uW6YKW+7YWQsq+/3Jtr30lll1X/\npOkLSt3nl2/bE5xm2qixvieSv6WKVELKzupmD7VydFVZ6LUl7oEC1iHbf7xcinpdqAqThaMf4WjJ\nnZ57usbqc3vzcMqjepM/e1Ln+HfIaHwpzoEvq7y367RT5SM3ZuyFhb4BBeyu3sOda5nXdRbh5hHY\nZUO62+kcoqFCNl0d9UzLdjsutODw5L8ztJB2sJ5OyCULtYIC1ifth7ZnUtgJ7Xrjt3CVUemMQu8R\nTzTz6NfsC77fb/r30ivbFx39JV1jjQu9sB5wqj35L6d6faGjkv9azkQn3krpV3pjezDXNsnT+LdE\nJd8EnTK8Fsn/WhyB2fbanhdA8oOeCwDAJG4lBQAwaashRIYUsLFSepP2eK19CljpC4/ABurf4SXt\n8U77FDCPWzDghUh7vNNWBazl/itAL/HCUDkyA+4Qz3y3WQHLUmsDtXEe4qmzWwnUViPx1KmFZCLz\nuQoRAGDSPkdgLfd0AYxqv6EX8B5aB63DqR2VwxyjKWQ0bOgwkUIMIQIATKKAAQBMooDNpnZtD/Fg\ngjnN+vn8+W95PF0EQ9JHAQOwA7///37//Xf8BRujgM2mdl6UeDDB848t9ksJl/j79/nxXCAYkj4K\nGIAdpPt/jsO2RwEDYFulRHFUszcK2Gxqp2qJBxM83ayVQvX9ZiqcYJoJhqSPAjab2kg38WCC5c0a\nVYfl8aQEQ9JHAQNg2OdzPk5IadgVBQzAKzBEtx8K2GxqI93Egwkeatb2uUYHYYJpJhiSPgoYAMO6\nhgepEZuhgM2mdqqWeKz4/ChNnRxPF4Vm/fs15/XxRARD0rfP88CAjYXPtkifcyFevR7yyg+NPzgC\nA2wz8dymh/R+7ux3wmDXiwrY6Q2qZ4UhEESAeHYVjjeGazUah5wzKfx3YRjpr2vD0JxUGalWs3nf\n7eicHs2x9cfFeCLHN6UhxGhHUxlm3E/LN8DGvvFVTKTQi86BMXqA/YhUrPkdRLZluFcNIQJ2fb/f\nY2znZ1BBZReePspkjjtLlFl5uIUCNpvOfscjHiu+P45fo6krgjqW/u+HUutJNavmwJjUKrLiXQWM\nUUTgOZqFoYRdwQY2LGDiHRm186LE8xL+KtzhG0fjBRFqzSoWjnN6q8iE3QqYePUCFpqwh5wzyMFl\nhPC2KmDZC7Gi73xk/84kJqWToheYFn6OXT7TXZxQ2ICBK/3b+QL2yX1j5uiyLe+7ZavsQsRTpxZP\noyjsMO3HbgLp3LLzH7sa73+En32FczLDiWqZphZP1j5HYJ/CF/6dWJdTLSeIB5eVCkm6xak1q49H\nKii1VWTCPl9krn+jk9zAa03uwD09NCfVH8Va+xyBAShJR/lMG9gf5UyYaRsWsPqR+PJ8VbsugHje\nZslohFqzqsXjJEPSt2EBAzDBTvvbnT7Lq1DAZlM7VUs8uKzUVuk4h1qzhvGIhKa2ikyggAE72+zY\n4rmPs9mKegkKGLC5tGe//EzwHU8cqHDwYxQFbDa1U7XEY0Xp5iAmbhqiFl42Hi7vMmef74FZoTbS\nTTwmZO8vczieECa79tQCy35VdG0FUVtFJnAEBti2ZMe369HCrp9rVxQwYAfp4Vc4rjj83sfuz0Ms\nM+/yBzTD77OcvnLg50q+7v3ITaL1J+mPSB90xxyG+PzcsjP8lGvv56s2zkM8dSLxVIYQfwpGXL2i\n21gnM7y1FbS8/e8dhAesxoFbbiWeVfsHkUw7qMWTxREYYF5lX1xhpJO9AGvGijcWsLVna9U6NcRj\ngn9OkHdcshH+Wxr2qXzXeCa1Zq3EsypStVVkAlchAjZkrzxkr/eQ5Q8ORIs3HoEBQMXx8FuIo4DN\npnZ5D/GgV+PRSThWr9asp/HMP/xSW0UmUMBmUxvzIR5MIHUJomuOZ2ZNIfMvoIABQAYDifr2L2Dk\nH96pJfPZOuo4KBK3fwFzYlmoNtJNPBurZ/7M7UKtWdvjmfatG7VVZMIrChgA3EFx0UQBm03tVC3x\nYIKbzTq8fnTFM+dkGJl/AQUMQIeuqwEHjr+t3b1/v8ftiVeGgchWBaz9PsrLn/0DwByuS1SzVQFz\nzn2/X3/XuNWBFKnFRjyYQK1ZL8fzXA1TW0UmbFXATAwiqwVJPJhArVnvxHPUsLEVR20VmbBVAfPS\nx9gcXZuoj6PzBDkmCU5qHI5W1bo3NPwR1/GnxNwDZQxdDDyyrN0nebLf5/Mp9Wy42zRapP0hNUeV\n/Zv5renduyE8/frsHNzqizgqjgImG+E1+pnv9jsC01/jap164jHN72W2P/U7fLMeuK6io7HLx2TK\nzSdrn+eBfYLn+zkLlQx4g5dc8Rvub0o1jH3ScPsUMCsVSy1O4tlVOKIeDgelI+3Zt2TfdUyszLDw\nrpMwRkUoO6nhyOwrVukNbIkGRjnv4BwYbhI/ExDtOrM/n82hY0O4ttXc3NbYVJcQz3xvt3NgAICX\noIDNpnaqlnhM89dueMr95cvN+lA6CKaZYEj69jkHZoXaXoZ4rFN72HHW/S8OjyWYZoIh6Xv1ERg9\nHmDObpNtDU94bwEbu93e/AoIsDEOLfCQ9xawgY47BbTc6FNtpJt40GinqwEF00wwJH0UsLui+9zw\nwAUAmIMCNkDULa33UtVO1RIPZD3XERRMM8GQ9FHAbiltYC+5fQ7wNPbqqKCA3VXZwKhhAPAcCth1\n9fpUKmwTTtWGl0SeXhupdupYLZ43mDNgINWwgmkmGJI+vsh8y+n4Rnrh1qMj3aUnJx01LJ2kNvKu\nFg+GUBtUF0wzwZD0cQT2oMkJWXnuX+Ml/kAJmQNBFLCLup54O0HLU2upYbhjfoeMYxLUvb2APb03\nz43mjV9k+zPX0xqmNvKuFg+GUGtWtXicZEj6Xl3A9ujftVcvj+MwAHt4dQG7rGvvH52+fuJUbe8s\nwxqmdupYLR4dx2NTGv8uRa1Z1eJxkiHpo4BdJJJsl88TcBxmi3/cl5fWqtLf2024SlDtQkRsgAI2\nyROb7s15UsP2QM8dr0UBmyHcw4wd6rm575pQw7Lfqq58vVp8KExW+kTmcFwxXKvReGPLpM/n38Su\nd4WTSiOf5Xe58qTrYaS/jprhTpP0R6QP0o8hv+/0AocLQ3DXRu2GXxM8cIa9l4F0zbYy59MXKEgL\nw5yFhr/64cEjjCgk/+K0evWGXcmoUcnW9eUT5ax4gyWZ32urO3Ecm3243kWawJ8AGFtyRvGxPRFe\n4/fSnqigpVUkkg91vTsO/R0N8IR9Clili7qlsZ9vVA27cFA1vIxVYgjHLW0lSHiNhs9tn+TRAN32\naQ+E9j8HFo13Z//eNSmaedcMP5/acHzjDP3/R419B/9+SjNvnqFzP7et6g0je7OrS6v3N4b0Xc59\nwgVVl3W8QKUkHFchHr+GfwwnXfbciY+jlUXOrQie4xEMSd8+RyrZI7DTQ7HeY46bxyh+N3tznY89\nPZCcUHGu/+hk4Bjg5/PxNaO3Xbz774r+aPRo/lrYpdQacWh+bJJNt8B2Dx8iCzarWkhq8WTtM4R4\nzdhTU43LvPPm4dFGOXqMJbq2PcjwEbljN9cYw50Asodi12a1sSEp139W7+4Sz+Yv18aCIel7ewGb\n7OYQypwxhnC3Xr+AM3z9czFklxKuiiHfJQBgzj4FLD3LvaVpnyw8FAuXO7BytMRwLDQt3vs2MoAm\n+xQw93zdGnQAdHFs+aHDr0o02TND7uHKkY2HWmXd0awrBu1r8egQDEnfVgVsgiEnAy6UoudOazdu\nMxOP/NiGl3muurQ365zyJphmgiHp2/8yekHXzoSR3phP4WgJKKGALdNew/h+CACkKGCtRlWR4xud\nXct9rhes9vVJtXhe6IkWUGtWtXicZEj6KGDONW+xw6vIquUCJaUrd8YuojRn9uHoQgGbXR6CWwE5\nd7bFTtie1U4dq8XzcmO/nz5ziWdLkUszwZD0UcBWqmfshBvqAFnZL94Baihgi0XfFz5QvbDKkXWk\nH8RRwJoMvJg4PVUbjSV23YrwiXjWUotHx+dHaerAZR337B8lCi97GmzmJfuCaSYYkj6+yDxbdqS7\ndN+mVfEspBaPiOzDFsKpK4LqoNasavE4yZD0UcBUkL245hM82dKW8JDLYPhYjyFEYE/heOPNZ6g+\nMck/QDUUFDOJCF87qTJSrWbz20emgy2FlzlXPgYaOzTfGNI0xFO3JJ7kjNG3NISYvjJ9jYJsPMd2\nN/+SJbX14/RCUosniyFE55SedA64vm9NiVasFr1PTwUiFLDZ1PYyxGNC+rg7WxWr5ZE9MwmuOsGQ\n9FHAABuiHVz9V+ANuIjjBEOLAKCJAnbu0W90Lkc8mECtWdXicZIh6aOA/TPtOg61oR7iwQRqzaoW\nj5MMSR8FrIYuEQDIooCdoFcEAJq2KmAmvkOuFh7xYAK1ZlWLx0mGpG+rAuac+36/d+4LF76PdAIA\nZVsVsJtnQdN3PzF+qHaqlngwgVqzqsXjJEPSt1UB86I7FISDirJ3z2SS4CT94Wjg5SzdjSYU7VmO\nm+u4vx2ZT//tdvyte2c+Ww/KLqSQAqNhQ4eJFLJ6BPb9K/z7/Zk/2u1W69QTDyZQa1a1eJxkSPr2\nuReib/4jCS5Xsqf7HGqdGuLBBGrNqhaPkwxJ3z4FjOYHgFfZp4ABeyuNLtwfdQCMsnoOzC61kW7i\nMcGfUU+/5nhcuHTn648TqMWmFo+TDEkfR2CAeen1t8AbcAQ2m9pehnisyx6ZObFvQIZfdFkYxvHr\nkWY6Xz10gbUrKnqBMgNX+t/xsfBVBihbkkLR7sMXp3Cf2/UzcIGJFGIIEZCjv+MAFDCEOJvasTnx\nmOAPwrzovjPR3zWpNataPE4yJH3SSX+f+FYNfUZTyGjY0GEihTgCAwCYRAEDAJhEAZtNbaSbeDCB\nWrOqxeMkQ9JHAQMAmEQBm03tvCjxYAK1ZlWLx0mGpI8CBgAwiQIGADCJAjab2qla4sEEas2qFo+T\nDEkfBWw2tZFu4sEEas2qFo+TDEkfBQwAYBIFDABgEgVsNrWRbuLBBGrNqhaPkwxJHwUMAGASBWw2\ntVO1xIMJ1JpVLR4nGZK+DQuY+JG4WnjEY0XlWe/6z4BXC08tHicZkr7dnshMEmBL4cOZogc1hY+4\npBePV9mqgPkNOKphgiVNLSTi2ZLaaiSeU4IhiduqgKXokOIN/I4vzHYyH29gtYBFXZXjwOv4lw0Y\ndqXpXXllaWgR2J7VApZuqGzG2AYJDLSwWsCAVwlP7oaXbKR/B96DgxUAgEn7H4GtHVGs9I7D8xyr\nIlTovJdiUFg/RyQW+3mama/QrKR9I/3M37mALb8m9fQEu85+eVWm1mNYvvEsT6Frloet3KykfYvl\nKdRowztxHPwZgtVR1OjfQGGt5etHP4WyxMNe3qziFNaPeAoddi5g+qKT8IiwfrZEs9axftptMoTY\n/r0ZnRhMdHAWYv20MJf5NGsd66fLJgVModW7YtA/O7oW66eRwlpqj4FmrWP99Np/faldixXeejWa\nNJ9gDFLrxzO6W5HKfKlm1YwhupWrQsrpZ756fAAAZHERBwDAJAoYAMAkChgAwCQKGADAJAoYAMAk\nChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoY\nAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADAJAoYAMAkChgAwCQKGADA\nJAoYAMCk/1YHYMbn8/l+v+nPLW/0P/i3RL+OCunyHI6fK+FlF9TySv+a7JzDRWdfk31BuJTTF5Q+\nNbpczvz09WE+jArp8hyOn8NcimabLqjlldHHPM38lk2jFMzL054C1uFyfpT24DezbUiyljaAbJHI\nvvd4V/Ti8O/+h/DXdL9W2VGGcy4VzvSVGOVarp7mw2VjMz+beJVESj9FlJ/pCyqZny49u2lEMaTr\ntvFTb4YhxA6nO/TSu6K/RFnb7vPj+DX6+4Xw6u+KNqdoUjSf+icqTT194+kRQLSneFsndIILmT/w\naCDN7WO3fjnt3dXMr6dfaValABozPw0m+str054jsAEqY1npC+6MxrhC76/UE6xHFf59SBe7HnPL\n67tehuUacyx6vULah/NsDya7iHrMjW+5FgkoYH2yXdHTDfLyppJqn0nLbmLISM7Rkw1HArP7mnQ9\nRAdPx5YfHlxG+6lwKe5q6UWvrsxvzId2Y9O+/WWV90af0SV1KzuEeMwk3CjStHckdhsK2ACVTt/Y\nLOw6lGk5LrwcW/Te031Ty7Lah2LYtkWUcmxsS41N+ztRpW+MPmkaYdcIOXpRwLqlaVofpjheHPXC\nerM2nVUpnsb5pzMsLSL0+TnPXH9lywdPq2D6FjZ+HV2ZX3lvV5MNT/vsPFtii0pUY9pn35IdeCi9\nABWspq2Q93gh0v61OALbweWjOsAu0h70XAAAJvE9MACASRQwAIBJ9s6BVQa+GRPHxkrpTdrjtYwV\nsMrXIUvf2Ac2UL8JBWmPdzJWwE5duFcNYB1pj3faqoClXVRuL4b7xAtD5cgMuEM8891mBSxLrQ3U\nxnmIp85uJVBbjcRTpxaSicznKkQAgEnGjsAq9w27fLM1QF/pRpGkPd5M66B1OLWjcphjNIWMhg0d\nJlKIIUQAgEkUsG19Ps7CWVhgMNL+PShgs024ticsXadlTO1aI7V4MMScZvXZ/v3aS3snGZI+Yxdx\n4JTfCsKx62OrBjYWZr7/l7TfHkdgs004Lxot4diYV8XTRS0eDPF0s6b9tvDv8+O5QDAkfRSwrZQ3\n17lxANOlSV7vumEDFLDdVGoVWzK2VElsatjeKGCzPXeqtj7j0pasdupYLR4M8XSzVvpt2UmCaSYY\nkj4K2GyPjnTX552dqjbyrhYPhniuWRt3+9HLBNNMMCR9FLBNtPfe6OdhM6d7fkrDrihg+2jZStmS\n8Vp03fZDAZtNYaQ7DEEhnpBaPBjioWZt/6ZX9DLBNBMMSR8FbAddmc9BGF6LGrEZewXs86Pygpnx\n9HroVG3vXI+VpHbqWC0eHfXMf2fad8YQ/rw+nohgSPqM3Uoqemh62uTim7EIf7M4GFLP/Hem/bUP\nzf2ldmLvCKzCxANsgLHenPa9n/ut62lbWxWwrHDUJeyoRqMx0yalL7s9Q+dc97uOgzD/l1VrQ3lS\nZbzOBKnMD/9dGMbPr79TpVJOZJKhzDfWd/uUB1KiNR4+dn1mhPNdHhJhLKWFSAqVMj+b9k4m7EeR\n+Y8ykULGzoFVKGy6wZURS5bfxx+EmQgVJQpp71ZkPqkLZ24I8fv9Hke4pU7oWtw8FE9IM18n7X0g\nhsoJFzFtw1gBc859f4R/iV4wPSjnmjuhY/c79/uhaluyzn5ZTZT5Imn/s/R//5Zaj2Y9xSq6wF4B\nUxY9B1mfoV4zMJaJLRR1byxgn8/juVt9uINc0ZDakgXXzx582g9v63QMILuIsc16/1MIZhmZf8Hr\nCtixvY3dkqVqQBe2mveY0NbT0mnIguxutvBeV8C8JzazdJ5Pbx5sfmgUpgpp49F128C7CtjMTbe0\neYw9VXt/I/x8PlIXZXEq+yHHdRYDtbeVWrMGX2ReG8gvtVVkwrsKmBO4zkJtpJt4cEfaXNn+kFqz\n+nikglJbRSa8roAdBmbLBj0nqYMwjLVZy/IVZhzeW8DGYouCsr9PEtmtpN3B2jDtRQVsfppmt41R\nI92j+qFHPCJbMmcCtqTWrGrxOMmQ9L2ogLkVFwoCu9pm2xHpuuGCdxWwyJJxP7VTtcmjEVcF8o/a\n+kFFqa3SkqDWrGrxOMmQ9L26gCHCFrSf5T2SsZ67gmOzFfUSFLC79rsmii15M41Xur/ZZpvwe9gr\nYJUHhgo+SzTdUwwJb+BHTJ6IOGzO10g1n5RSegumfUotvDSe5UVdbRWZYOyBlpUnMnulBzHrHCeN\nGuke9XFy63Dl6uJMQFY98/WfP64WmFo8TjIkffaOwCquZQD9Hpi2ZMe35Vaz/CAMvbYqYF7UDw1H\nV8KD9J+H27rSpPBXtUnpK8cuK5wo8pHnT9IflwvlRh1qmR/+emGSC8YAsu/yxcBg5v/WMKlsnDnJ\nUObrjjlkfaoDKT9bzjd9fWVM7OZw2enboxekYT+x0J5Z5ePxCTw/O4asn4FE4qlkfpr2riHzn077\n6DVW0n7sUrqIZNpBLZ6s3Y7AJq9xnVNrw+36ubZ0YV/McFkFa8YKYxdxfL/f49g26pD6v6dTl4uu\niZDqh7qzeOZXaJ2Gk5JmvnjaR9QCq8SzqrSrrSITjBUwl79q7pv9O26iky4le+Uhaf+QjQdXdrLb\nEGIWuXgHNQyR7VPC7y62/5gbeEUBq5t8nKF2eU89nvlbstr6QVZLjzDcsu4369hu6Gk8Sy5fmr1I\n+yhgs6mN+ZzGM7mGqa0fDKHWrC3xTO7aqq0iEyhg1y3pMC06vbxgoYACjouUUcBuad+zD9wMVpUT\ntmRbWtqLNq3jZJg4CtgM0ReZ1wWS0RjPtC1Zbf2YVu/rzOwJqTVrezxkvjIKGFrRG8U7MYQua/8C\nprbDVTtV2xXPhBqmtn4wxM1mHf5NmN54fm7tODKGZBFkfrf9C5ijAzUUx2Ev11VLdvouPJkv6BUF\n7JSVzUwkSLZkvBOZr4YCNtvNU7XDjyavxfP9PjWowqnsLak16+V4nqthaqvIBArYRRcG5X8eUKI1\noHknnmNjHvolAa31gyHUmvVm2j/Re1NbRSZQwCbZNTn9xuxGlzHc1ppwtNo1pL0Ce3ejr9B/qMSu\noo2Z1T9TKe17b1coS/Zm3GkN04xzY/sUsPrDmnWoxTYwnuDxu/m/T47nDUj7ax5K+zDz+08xaK0i\nE/YpYAUkxAJ/7zxS6+Ozwdpy4XDNxEHeEKT9fNsXMOeHWMKH2P7+9e+oyzFSEU2qvKtrkq+m/im6\nl2YY9vWGRXismSc+8s/s85MKG/lXbH9ndWcTtk4987Nvyb7rmNib+adhjIpQf1I5w8n8bvsctGbH\nUrqOyruG2q9dhXhnZcueCdib+MBOaQixPeyn0/7yu0a9HdeIZ77HVYhTifWwAMCwfQqYH5rznu44\nXOsShuOTOojHtJlpf8flZn3o8EswzQRD0rfVOTDlDfigFiTxWHdzjfmLLJ5e62rNqhaPkwxJ3z5H\nYAA0vedCRExGAbOB89hYjiIENRSw2dRGuokH7bbpRQmmmWBI+ihgf5BCAGAFBezXrN6lVidW7dSx\nWjxY6LmRc8E0EwxJHwWs252NihQFgFEoYAAYPIdJFDADHr0EUe3UsVo8bzBlYECrYbWicc5JhqSP\nAvZ2aiPvavFgCLVmVYvHSYakjwIG4ByHBxBEAQPQZPIRAl/exykK2B9T7nmjNdatFY1ePBhCrVnV\n4nGSIemjgPWZ3yukHwoAWfbuRn/0U9JznpVJOvyTiHWorSu1eHSU0ntU2j/aVSLtTwmGpM9YATt9\n5uyFBzED+uqZfz/tJwyez3luC15lqyFEihZeaELaSx08AYetCpgX9UP9w2qPn7N/b5+UvrJ3hunc\nbs7w5qTj37VhCE6KXiAuPfwanfm3VnK6Js/e9Xus9kRbh2tGJ+VEJhnKfOmhtmgl+qenVwZS/Ouj\n6tX7Af0yS28aMgZSX8QTS8RlF1JoyELDX+uZn6Z9+pq2hT6b9l3zIe2XW5L5vaTPgV1YfffXeOVk\nwKhOSfv5BiPdIAzWm8b6OxrgCdIFLOW7osfP/gffU4hGKh7apOfvKNg1weUyf2baA5qMFTCX20T9\nX6xsumoH5sRjRfbKw4Hr6tGrBH2z6lyIKJhmgiHp2/AiDgADzR/HZuQcjShgs6l1sogHp+63Sf9Z\nvbtLPJu/XJoJhqSPApY350uda2MAANMoYBnZntCSsXv6ZABQQgGbTe0bgsSD0EOrP7hOUmJoQTDN\nBEPSRwGbTW2km3hwSNf9uO8+tjbrnKEOwTQTDEkfBUwUvTGIYL8KWRSworCEPFFOTsdS2HEAQAUF\nLC8tHqPKidpIN/Eg8kQLqDWrWjxOMiR9FDAAv8KO2kN71PrtRhl4QDsKWI3fzMZuxi2namd2xdRO\nHavF83KjWkOtWdXicZIh6aOAFfl06nr0yYVFlGoVyYxVfFoyoAV9FLAaX0WoJXinRzOfAon77BWw\n0weGjj0XOnwbVjtVSzxW1DN/eNqPzfzkEZ3Z18zrLAqmmWBI+ow9TqX+RGZnIQlyMcdnzmce86mN\nvKvFI6LlWeTK1JpVLR4nGZI+e0dgFRYfqGMtXsixmPbe01+1xPa2KmBZ4ahL2FGNRmMWTnLB1vt3\nk1aJ8J2TKuN1JohnvnPxuj2qsEiEr51kKPOl+27RSvQPUC8NpKQvTl+jIBuSHzZ89IrHrngWIh7X\nk/nZtHdGVuOR8PMzX239OL2Q1OLJkj4H1rX6ZDfdRkuqF9B6CQ0AAARVSURBVDS1J7DptA+v17cW\nOyRIF7CU74oeP/sfbG262VAXhq+26tTiEZFm/gZp79ZlvuCqEwxJn7EC5nLNHP2FPMCW6nlO2uOF\n9r+IAwCwJQrYbGqX9xAPJlBrVrV4nGRI+ihgs6kN9RAPJlBrVrV4nGRI+ihgAACTKGAAAJMoYLOp\njXQTDyZQa1a1eJxkSPooYAAAkyhgs6mdqiUeTKDWrGrxOMmQ9FHAAAAmUcAAACZRwGZTO1VLPJhA\nrVnV4nGSIemjgM2mNtJNPJhArVnV4nGSIemjgAEATLJ3N/r0cSotkwDrSulN2uO1jB2B+WcgedGQ\nsf81O0mKWmzEY0Ip80n7a9TicZIh6bN3BFZ3bM+rAwHmIe3xTsaOwOpKXdTj19LfZ06qhLckwuPx\nvmvDSH9eG8bByauMSehk/pFma8M4fg2f507mR5OsZL5zTvqp5NFK9JtomHlh8NlJH1OPXYegJSnU\nnvktPwMXmEgh6SFE/dUHPIHMB1oYG0L0XVEv7GyWJglSOzYnHhPS9Cbt71CLx0mGpE866e8T36qh\nz2gKGQ0bOkykkLEjMAAAPAoYAMAkCthsaiPdxIMJ1JpVLR4nGZI+ChgAwCQK2Gxq50WJBxOoNata\nPE4yJH0UMACASRQwAIBJFLDZ1E7VEg8mUGtWtXicZEj6KGCzqY10Ew8mUGtWtXicZEj6KGAAAJMo\nYAAAkyhgs6mNdBMPJlBrVrV4nGRI+ihgAACTpJ8HlnX0U7LnPI9nTEyNqYdabMRjRSXzSfteavE4\nyZD0GTsC83f4zz5APZyqfDCuFhvxmFDJfNL+ArV4nGRI+owVMAAAPHtDiHXpWIpgv0YtJOKxLjuE\nqLYaieeUYEjipAtY1Jz1MeLw+aHHzwwrw6L2zM+mff0twDakCxgbId6JzAdafMxtKum1WEfHs36B\nImBalN6kPWCvgAEA4LgKEQBglPQ5sCHCM9tLlu5/KH35tDR1DoXRp1IMCuvniMTiQIVm5is0K2nf\nSD/zdy5gy69JLV0hdtDZL6/K1HoMyzee5Sl0zfKwlZuVtG+xPIUa7TyE6G9PsDqKms/nYyVRlli+\nfvRTKEs87OXNKk5h/Yin0GHnAqZP/w5Aa7F+tkSz1rF+2m0yhNj1lWeRGEx0cBZi/bQwl/k0ax3r\np8smBUyh1bti0D87uhbrp5HCWmqPgWatY/302n99qV2LJfX9U8EYpNaPZ3S3IpX5Us2qGYNfRQqx\nHfQzXz0+AACyuIgDAGASBQwAYBIFDABgEgUMAGASBQwAYBIFDABgEgUMAGASBQwAYBIFDABgEgUM\nAGASBQwAYBIFDABgEgUMAGASBQwAYBIFDABgEgUMAGASBQwAYBIFDABgEgUMAGASBQwAYBIFDABg\nEgUMAGASBQwAYBIFDABgEgUMAGASBQwAYBIFDABgEgUMAGASBQwAYBIFDABgEgUMAGASBQwAYBIF\nDABgEgUMAGDS/wHK2/e3X/O9zwAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "% p23.m - eigenvalues of perturbed Laplacian on [-1,1]x[-1,1]\n",
      "%         (compare p16.m)\n",
      "\n",
      "% Set up tensor product Laplacian and compute 4 eigenmodes:\n",
      "  N = 16; [D,x] = cheb(N); y = x;\n",
      "  [xx,yy] = meshgrid(x(2:N),y(2:N)); xx = xx(:); yy = yy(:);\n",
      "  D2 = D^2; D2 = D2(2:N,2:N); I = eye(N-1);\n",
      "  L = -kron(I,D2) - kron(D2,I);                %  Laplacian\n",
      "  L = L + diag(exp(20*(yy-xx-1)));             %  + perturbation\n",
      "  [V,D] = eig(L); D = diag(D);\n",
      "  [D,ii] = sort(D); ii = ii(1:4); V = V(:,ii);\n",
      "\n",
      "% Reshape them to 2D grid, interpolate to finer grid, and plot:\n",
      "  [xx,yy] = meshgrid(x,y);\n",
      "  fine = -1:.02:1; [xxx,yyy] = meshgrid(fine,fine);\n",
      "  uu = zeros(N+1,N+1);\n",
      "  [ay,ax] = meshgrid([.56 .04],[.1 .5]); clf\n",
      "  for i = 1:4\n",
      "    uu(2:N,2:N) = reshape(V(:,i),N-1,N-1);\n",
      "    uu = uu/norm(uu(:),inf);\n",
      "    uuu = interp2(xx,yy,uu,xxx,yyy,'spline');  %%% changed from 'cubic' to 'spline' to avoid matlab error\n",
      "    subplot('position',[ax(i) ay(i) .38 .38])\n",
      "    contour(fine,fine,uuu,-.9:.2:.9)\n",
      "    colormap(1e-6*[1 1 1]); axis square\n",
      "    title(['eig = ' num2str(D(i)/(pi^2/4),'%18.12f') '\\pi^2/4'])\n",
      "  end\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAgAElEQVR4nO2d27Lk\nIK5tXR37/395nYeM9nH7IoSQQMAYDxWrMm3AtmAiIZz//v7+DgAAgNn4z+gGAAAAWEDAAABgShAw\nAACYEgQMAACmBAEDAIAp+b/RDQA4/v379/uDnFiYAiw2CXhgMJh///79/f0xEMAsYLF5QMBgMAwE\nMBdYbB4IIcJ4fgEZxgWYBSw2CQhYFL84g+3E3x9fp99Kvv739VzhgOfxXwcIp9xq0Z9yRWgww0Qf\nWm7186HfPjdY4O2YrwOKVQgl6Fv1RLB/LLYbCFgULUb81ROuXffrv7dzb8c8D7j+97WE54Uoq9Bc\n0Xn665T2WTLE0Tjsvj7Z14mIYIHPol6F4VbXzYZvlQolfJ3y2qprkw4sNgcImBvPSdx1enioB4ji\npE9wko43GRPa8FrX63T1a+x4bYZ8yq098kyWESGOV4s9Hg+oqrTi8QYLlA8w6G5tI4/SPbkKIRbb\nEwTMhy/v5zXu8fyvPHvVo4y8CQ1TTkXl8oXLP28UkZaxvFrs8+/b8SdVTvZr7V8W+CzqqWfFGMBr\n419jm0ILn6djtNlAwDzRT76Us8v2xpz/FscCZTuLzXv9qtapem08uONiscVvnzUKFvg1T9LEtI/L\nqtjvj2IYU76iW5nF68JiO0MavRt//0Vz8L//5VmUb3ue48VROQlVNu95ilzL0zE9G/zaeHDEy2Jr\ng2aCBdba5GsJ8kWZzV5ow9+FA4vtCB6YD39v2VDCxO3LxK/zuON/HRdbG+TCb8c/P3n9r1CLpsxb\ne2AIX8/ly2hln/vVYr/MqcUCIw4onnIFo00Frm4niCrAdGC0kBw8sFhklwggIRgtzAIzLAAAmBKS\nOAAAYEoQMAAAmBIEDAAApmQ7ASMLFuYCiwX4YqMsRAYCmAssFkBmIwF7vjyGAWI35sq5fX3dEUa7\nFXNZbH82ErBXBtrHwF2i21Y9pF539nx2e1Y9pN6J2G4NDAAA1gABG8ZA52/PqqGRPc0Gi83M7iFE\nAJiFtCE1RG4U2wkYpgZzsZvFyiqV827Uvr8fvNhOwODwnsnSUaGFmzXOaE7CrxcJB0A7CNiCaPTJ\nsUcp5ZA+DFeWH9+vF8UL/oNAwBZh4HCg/0lf24mwGD9L2Orp336EeqtrDwUBm5uJZnbFH7fNfwnQ\nyETmGsR1c/q2N8ERBGxWFugDt8YvH1PanJ47gh1XeSPafMoYdt4IAjYfC0jXK6wZrEqoxUZnAMaV\n//f3h4Y1goDNxKrS9YQ1g2VwH6M7R55DMwzRsEYQsDnYdhxnzWBqvEbnbOFlogVJQMCyQ/c4/jtR\nZa46Fy7PK7/9n3Ms2zQLJ6wFBCw1eB4nuGJz0Tgo59etGy1BbzTMDAKWFEbqV0jfmoKWBzS75TPT\n6gkClhGsX4YZa2Ya1av/O2K+aMzOOJhpxYOA5aKbdLX07Qx9Eg3LifmhuFi+79te2tMdf1aqORd7\ntoGAJSJOvUI7dnuBNujz2WhRLy+DdLSHr432hvUtrybBDQQsC+7qFZd5XNwZ416jABo2O41O29HL\n2MxpGsy04kDAUuBo36PSt247Y/ps3GF6mwebDZvP+v0xRBWuaRrKNmg0DJ0zgICNx2UITpV2/Cpm\ncQ2j20+K4cHlyW+6bv/y0jCoBQEbTGOHTKVbr9wCL0dAPJM9zsOxSVHVKXmk60pt0jyG6gsCNhKX\ntetZ+kPc6w0JJC5PclNXJs1jqO4gYMPYdrMnOz0XI9r9avdaamXDVh0LXf1BwMZgnogtM+777vRk\nXJiIKuNvn+cd9f3FnIKEHXYGARtAS+bVYn1Dv9MTVkL5uG3zvPbV1lsWUlVpLRqG/tWCgPUG9brh\n6IrR+ftjyMWoUq/hiR6GtVtZh1ApRxCwrkydNxxKe69mhXw9qtbJqo43NwbtSQUC1o+c6lUc9Lt1\nV2amcNI5y0OPMuiNE9YHBKwTHTbKKMu8fSJXcXunhvIsM3TstVE+3LTq9aMqaR5jDgUBC2d4KL/l\nrU5fx4e+aJE1cFAy9nFjbMNBwGIZ6HiFvqSjJU1LUzhDAxQZayQaP0y2ZIy8HQQsEPOiV2Olvz96\nrgpcq2ZfF7SgDzN2aEwRs6GSc+QCApaI2d8s5f6KDTQsPwOnGlPYBlOxUBCwKGoXoo/mV/om6SRe\nMmabojJe7IDL6zm+MKwTY3KjQMBCMBi0eXez7dxoXGSMoQEcUfYXQ7fCUEeBgPnTx/dK5XV9wd7P\nTej5iKP3U55GO+pFnfQXPasJmDB7isv8vtUSrV5TSNcV5d5P4dyJLtbAl9H2sdh2kucjmJMsDr/X\nR+lPSX4zs7GUgF2t4dWYokcB1OuLv4afnVxbw2SjXfWqu9EoBi12Cx34z+gGdOX5XomvN03YCg9V\nr7Opk/alv7+/UbNLx6fcmdeWJ7wc5ZPtbAC/LtnYX/R2Kx9WdeEJH3FOlvLAijyXZLzEIDrGPbV0\nXbG5U41O2DUcNBevi4gLmMFu1Cr3vBbbmY08sFTrXobylxm2oqfhy6wizPXEs91237hftquDH7t4\nYHFR7A6Rw4iWK3tj0E1be1nLhVXvz7yPXtPyea9uUpYSsOss6bYw/vpVO7adtlXlB7VWWeztAn17\nZtXVrTo0PC0z1GJDWfUZQVqWErDjrau7L3edVLlTtb6X16JXYyr2V2J3e8PI7zr5yjyc8c4oNUw+\nJq1VoNDZWE3AumHLRIgr/LWQ2nqLXP3a9pINK9sMHwsgP3dWm0DPRkkcA4nOUXyWcIah4haxTheq\nsRxGq8VY+5kWr+7rgLVvyyjwwCzECVKjiXfOtn9N8oaJKNqb7ckWfWWcaXABAasmTpBa5GfgRrFz\nahk02IE7yvDyqL3zskmsZDArXcsQELA6Qn2v6aTrpN0VU56r7/CMC69UWUvLkidOGHSANbAK4kLY\nLeoVt9BVi3lVzL39SW5INszWcl3y1D/fllWf9hUj1px2AAGro3Mm4fAqasnWHvjhsmPhFD8vDTNr\nDOIEPxAwLamCh5l3UNkGl6qzGLyq8HXTaxNQzRpWNAlNBgqmsjYImApDLob+YIN6Hbl9ndCBI/OF\nJyToQdS6Yu7NwAzgQMD01K5gaw6zqVeeRS8Bg4YxX3Yn1Fr0GiY3ACcMzCBgZYI8KrPvNQuMHWPp\nEGTWhxNbFsMEeVO2E1YFAVOh16S4Nszie12p1TA0z4ueS6SnjBUPc89Qjav3i7Rrz3uCgBWotdcg\n94tucwWdkxliLY1a0hJIdGkbzAgC5kbVgnateplalALGjs4MvNvKZ22IFrZXaisckoOASUS4X7W9\naMbI4ZM+u19huLUUl8SEeKP8VeN+MmWcU0YzIBAs6QmvkvKhymqrtohZW2Qv2b37oUlbcWYn+hrS\nn8erp1oKaaz99XTUrhEE7JMI2zKU6diGq4p8FXudPvtePn01mlR3WJYK4Vv3r16PPEhiXAIEzAFl\nz7EFD62Neqla2b1vZylPLBaLExZKe3BMc1htRF0wPOFbWY1cNOyo6RQIXloQsHeqtnPpi60KHrp0\nmJa+d/bzzlN7l2DRhpgjY/pza+c05rWrLxvQFFjVtuLlV5kidtsZBMyBiEijVyEuzpNLab6yhMi1\nY3im58FVPn3RD6v6yhycfD3490fPJWFwBAHrRNB+Mq/qipzTVfpzQmwbMxqnI/pyZOMxzEUcNew8\nRX8w5IE0+heq4ocZEj3cS3jlOmaZS2AlbDiOKzpnbromi/34Np4vwyjGGIUasbQdQMB60M1xiXaS\nzgErqHwwYFikcffONTImG4+vhg0xVEPiCTSCgOWixdBXypVCJiOIntwcba+Hd/fDsKLlQcDs6LPn\nu4lKn4paxgXGlFF0uO16V2w9DcPHGgIClohG96t/pnu36qCRiMjhFxoZk+0nQsNCzRX1GgUCdmdG\nW0RLnmwisWmvsRhRFJytrxMNsqdsTChf9aZ9dhOBgBlxjx+ahXPU0tcmCpGf4qMfNScrZlIIGlYb\nS9TUFWSxmtsrZPy7t2cr2Ae2AnQD+KJxYvRKnz1bX18ZTrkdUHsJAu6TgxkjQANBwFJgnhiONffi\neOF4FvThNEXhAV3NVfMcZdkwWEKjhh1OcQuNGROlCAUBC6RqmGZAB3dqt4gdakE6T1GedcpG1U4p\n4awWDTs8ZEx/b+nacSBgFpLMqpI0AxbA5hNflezQyZisOl+FfJ3VomHH/8qYpv1npfqDIRQEzEgS\n883QjLh4IJFGGfn+6POMjmZD0js0glCd+e63r4r5igZFfJZziKHR2qjpeVaVxwm1bCRgtfOsbthM\nOUkHiMtFJMuxg8X6WpFSxr6E6hC1SpY9QcNssb5nMzJ0N7ixi4BdjbjP0J9EYGBSGi1Wo/1BJqqR\nsaKG6bVNOOUoeWnFq2iBEaADu+8DO1ehYWEWe8rKy5FHz+gbciqHoDqyv/XqAwkFynn/nW2gsa7F\nLDaO3QXs14vO/zJpWpLbU56d9sv51+XNUmcVRWl5/Vwo81Xb5JhzsSXutCxPLmaxcewuYJOC0MIs\nFKXlEF8cVeVvFTWsz8Iq3bMbCFg1vtaJrYM7LtmJvsgrWF5+mFzR9YA4GaNH92SXJI6rWWNekJ9U\nFitoj74QOa9dzt2oSkaXKzr+N5zoe2+VuTOONW7OLgJ2JBgFAKqIsNgq/6CooLW7o+ScwC/hkT//\n0jChotsxysYX0RfFWOTFRgLWE8II7XAPB6Ici28JUJpTDt37OF4/f56lUcRiusQZUWyxNxdzxeZr\nYQ0MYClqA25fh9WOpOe6lDLFX97U9fq5Y2m3Y6oaf+PLC/w6GIlyBA8MmqBDLkbjA9XH5Yp+mMt6\nWPGsr8ZfP/mqrngMRIOAAUxJra+gGcS9piNKGXPUMGE97Lj4YVULVLJDhm5lAAGbFVwfcM+g8y1Q\nmUnxJS2yVj1PKa6HFRvzeooj9Fl3WAObkqm7Ad04IXEPpbgKZVjcesb6lNVplsSCIHs+AgQM7CBF\n2fh6ImOfVDE/oltOh6YxcdBZ3EHAAHang7wJPtN5wJeuyLn1QmlCXZ1dMaZ6QSBgszIwGPIjeUgk\nefPy0O1GFWXDFks0J9B3c8WUVWCxBhCwKJTmOFyHWkg+qUzevDz0vFGybNhiifIpGlcsrg9WbZHG\nYmtBwELoZojTrUhX7aW1VQE3st1Mzf6wqq/kUzTtCZ1Hprr5i4GATczYjkG3XACbtv17UFtCS66g\neyzxepivjGWbOqwH+8DmRv+WAUfoltvyGhC7jvtVu6wEQypuAqsqULmR+ZZp0mLkVZFDMIMHloLG\neV/nfKoh58JYfir1i7bdvvr7L8f/ilmRom9k2OklxxKVrpjhWq4NMKgXXcMGHth4GuPvhlcMmGmv\nxauReIHt6O+h8sjzGP0IXnSbBD/s9XOhtVU95XYth+JyGr0u7NkAAmYh4ejZIZaY8KqhA4bnLmhP\nVRUtAcPa4OQXTyUrHgndIIRYjT4jtsqvcokhhGYDN/ZP9/k+3Ih4+uZnURWyk79tzM7wOvGvhKEl\n0AgClgIX6//1oohMKhf18moPCDwf08DZwGmQmiMFuzUsehU7QkRPMcOMzQwCthr6mW8RYd3eAF10\nOlwGVr2G2Qo3F+jYU2AUCNidUQbdmMpxK6pxgmnLpHKB2Wgfet7nKj+sz1fnAY79DvpDEsf/4GvN\nHRIr5NoP9c/LVh1ZC5oE7X3BlrWhqXdsP6V3tICAJSKiI11LK85GHetVVmoukz4fzRC/JKgLJNcw\nMIOA5SK0I431BTUwiKSiv5D8cPe09Bp29O0mGHwjrIFZ0EcaN4+w0z/hiqY7KGWmtnBl1Zt32OlA\nwNKxTBdCvcDG2C7QLcOeDtIOApaRBTTM0DnpzzvQbtuNTpgyH3KBPrgDCFg45p4wb/8JbTk6BzKy\nhunLCdUwzNgFBMxO6DA97xzQthJOf4Yr5gUtzf5lfbcKCidi7V4gYEb6TOWm07BR25+fzRjeBvgi\n2qo1+5erNMxXxjBORxCw7EykYeb3TtGlXTD4Fl6lDaFn87xkDFP3BQHrRMtwkH8oOeb5QbLdmMJ4\nXilGEc3nag74OqtFxrBed9jI3EQ3ixyyy1JJY8Po1TAE80sD/v73JW3KrdnKI6EKBMxO7SSu8S0b\n126TpyegXguT/AVLXm/fsF3gedZtELjJ2+1g8GU1AROmRcqX1YZ22vbCrxGMsb2C9+548WW0SouF\nsXw9NR5ZB5YSsOsA9zrYuZuUQZBcNOz3xxAZc+yfqwpSFbLRht6faB/li86B984XiEn3ZCkBK5LB\ncXHEEItvwffu1Y4ae6rdYhabk+SRUhDYS8DOEX+sE+ObkXFeVISSBUWxOufFTZqGd7xZ7GG6nMbF\nV8fxfd5n0RPukpKJBax2mfT1gPZuaeve7snNwgBX6+iYz9VXETpj+ArEZRgUqm6vkJ7wLEouR3hv\nheZx+GpYz4DkkAtsJ4/FJmdiAZs9ABXXZ4RcgNpz3bFdcsLHZ0N/FdkuuT1sEBHtcJ8IwlxMLGBP\nrtZ8Wxh//cq3XkOxvrFEoZa4wqtgrHnytMwWiw3NKnqNZyoJsnNHi8rmhIGGpQTseOsh5ydJRoHX\n5i3fc1rGr+VvzlfmoW9qnNdtrFVW88SxSla9WN7YFmM1AZuXhXsO6pWcWnl75g0d//t8J93BRkBy\nOhAwN1omudc124k6vAYUKA7fe2uIOtxO//qqQ2MEajsmFjsRCFgiFpOx9gsxDCX7jD617oI8jrc7\nH2vcdpywueBt9J64WP9rfGYufo3/ZSK0FOLYJCiSZ+zeZxYCjSBg75g7s5eGnXloScYUJS7SdSzk\nhvZHsMCicWbQML16xe1CGX4TQAkhRH98M76OGUbziC0Kma93YcZmk+N7QRUIWAi+o8BVxq6fDCdo\na5351jH8FdFYZp/tiU9qUy0iioW5QMCicJ/JnkVd44r9e2a0jjLctNNue523J9r0Mq5tbGqeBQTs\nE5dRIKIbfG24OcIUJbqKa0WMGtHozbKDK2YrH3WBHwhYLNFDgLAXR3O85qzMc/BnCXBSTJrXv+ci\nwuOPfkMH7AACFk7PaEyxitdRfvhw4KXxwy8kD8qEQ+Udu4avb59U0R5/ru1HqN3aIGASvvmEGUJk\n2TpzhnsyL0X7lL81mPdTyaoOaN9cganAFQSsH50XxvPjeCu4q080+5nMU7SvWHRcjo8hRIxVLA8C\nViAimXBgDmEGNr/8zvRMmg96oOZEj4jGQCoQMBVxCfHHTuN40PUy0f5C+VKJv/++8yXbbTQbzG49\na1sQsDJxr5aJS/FKxfIXmBllCCGbKba3pDGjdfgdAA0ImIrQjY1eKV7Z6PPqkM3HGv2bNaqS5kd5\nMO1dAN9rKxCwRLxmcE3XFXs2fnP10jMkaV6Pl81gD7uBgGnp+XaZ6wjS7UUYZvK3EH54Jc23P+Jn\nTL49w77l9KCiIBoErIKeGnbWePtk7E7kVPugmW6fVK1ytSfNuwzx7i/1iEizguQgYHUM34/8Kmnd\n5ox5OjbqdVKVZORiwHnuvHtnxK7mAgGrJlvOcZJmdINV+kauCYfz3saI9hM8nA4EzMKpYcfMQ8CM\n5Jk3pML8Uqi5El9D84MwrRlBwIxcE46x+w4wXYhgeN68hg55rfhek4KANbFGNCY53F4NLRlGcdmG\nNjo3gznovCBgDswYjckP97OW9izZYrZhnzWnzlm1GNi8IGBu3KIx1w9BD7rVgu9OD+UuDsfy+4CN\nLQMC5sy1SyBmSrhLjoTu9Jj96eByLQYCFkjngMwscB+iYWn2Bi7XqiBg/agNyEzd2YRLm/q6JmLn\npdkk2SgQDQI2Erlfta83NL7Se0jV4Msab4iWwaffFgQsL+2vN23RIYaAxbitzsoHZKNoyZkbD3Eg\nYMtCl4Yvxr5R0wb2DE+2EzA2LcJc9LFYOgXMyH9GN6Af2eaYAxuzZ9XTgcVuXjUU2UjA/v7+mGbC\nRGCxADLbhRBvjJ1e7TmvZErbyJ7Pbs+qQWZNAVPuAmF6C0nQ71vCaAFO1hQwOjnMBRYLYGCjNTAA\nAFgJcsoBAGBK1gwh6um2LezrlXSdX/DT+c14QnX932y0zBZAjHZIdUPexbWM0Qaxr4D1zCy6WuHT\nIod0hg4do1hd598tXACMdqzR9tSSZYw2lH3XwPJsssm2X7Ub3S48z7NuJM+F7Gm0Pa86z7POzL4C\nloefpW44HGx74Quw57Pb86ozs0UIsfOPA1VVt+0ka9sLV4LRZmPPq07OFgLW2fL01W27QrvthevB\naFOx51XnZ/enMjCh61d1kgyrbtWdN7z/LwUvMwaNMtpRz25sdQMt9ljIaIPg7gAAwJSQxAEAAFOC\ngAEAwJQgYAAAMCUIGAAATAkCBgAAU4KAAQDAlCBgAAAwJQgYAABMCQIGAABTgoABAMCUIGAAADAl\nCBgAAEwJAgYAAFOCgAEAwJQgYAAAMCUIGAAATAkCBgAAU4KAAQDAlCBgAAAwJQgYAABMCQIGAABT\ngoABAMCUIGAAADAlCBgAAEwJAgYAAFOCgAEAwJQgYAAAMCUIGAAATMn/jW4AwPHv37/fH39/f2Nb\nAqABi00CHhgM5t+/f39/fwwEMAtYbB4QMBgMAwHMBRabB0KIMJ5fQIZxAWYBi00CAhbFL85gO/H8\n+1bC61e3iq5d6xmp//rkWdFXCYZKX0+5cSvt6+oglJa7/WVIvz+qDMxgTs+W36zoLOTaib5q0Sxx\nCVeExXYDAYuixYirus3tv9ceeJOBpzA8e921wFchudV4vHXdZ09+nnW7BHn4gz60qNfXQy8aw9MC\nBRt+LeHVJp9tOx7qctre13zuVY2w2DwgYG48J3G36Z5+dBACFM9BwbfPfKlU7QFXbu283RN5JsuI\nEMerO242WrmWryoM5T/N6dZmpeVUuUryPblOB7HYniBgPrxONvVzxluXEGIRxTDF2ZJXwZBbbuY5\nSBWPJNIyFsE9en2OssUeH4akDAV/RaR9xUAjacLnGG1CEDBP9P1NHyT8+uoaxxAihMf3QOCy1PRa\nyNVLe22nptjbWYa2QREXi9Ub0pcxaOKEryXcPjwLvH7y2hGuFRV7itKfO7DY7iBgbmj8j5Ov+azs\nWl2/Kp4ypBc9Ky0uJ9xU+SqBX2eBCy4Wa67xqh+NJRyi5RRNqNbGvqq2lQYtIGA+vIb1hYnbl4m/\nlnObSH6d/jzg9t/iAZrG1FYqwGrBQL4eruCF6ItSWuzNkXo9RW9Otc2TD/iqFKNNBZOFTjAvg+nA\naCE5eGCxtM8fATqD0cIsMMMCAIAp4V2IAAAwJQgYAABMyXYCRhIRzAUWC/DFRkkcDAQwF1gsgMxG\nHthvl8noVgBowWIBZDbywJ4ww92NBfQAo92KBSw2lK0F7BhnH2O3iA6sfWzVQ+p1Z89nt2fVQ+qd\niI1CiAAAsBK7e2CjGBsZGFg7IZF5WcNsDG5NiyeEwYeynYBhTzAXWGwLtS8m7tAAHqgj2wnYkvSJ\nldPxYBbyvM7x1oB///7ladsCIGAT07kn1P40DEB/br+QmY3bL4cdiZs6BQjYfIyawQm/tKk5DCCU\n6fTg9ztk/GZNCwjYTOTsos8gyevnAHG4yED/ZA3lb0bDFwjYHOSUrldef6h+ipbDpJgFwDF40FLU\nzxWjjxhAwLIzkXTduCrZvFcBjRTdmkarMAz9EUGCxjgEGmYDAUvNGjZ988kWuCIoove/W6yiqoP0\nDG5j831AwJKypN2fEf9juUuDE4PzccRbxSir018dTpgBBCwja9sxMrYqLc/UkM6gPLjR0ly2QpOs\nEQQClouIzqakf1I+OcTL4DUdcfdCahumlCvb+zXwsdxBwBLRnkzllUbVp48xLV2ACE/ayyT05djC\nnoYSZA1D4WpBwFJgGwV8F6Wf7whwLFyul047KREP7mcP7eX0CTCe6LM2MHhHELAsVBl0dD5V/wx4\nuvRchFqFZoiXRU4jgUGXcMbGhcLRMC8QsMHYYvTdTP+5KzmiamKJc9HnSRVrMX8b3YmuvQYNCwUB\nG4lhF8soi9f0SdiBPgbQEkiUW9jTgFGpaBCwYej7Z6qM83NkCYq9JLlMeKXzA3KvrkoUXd4hIlv1\n8yt6QRUI2Bg6rzD7wkauDen/uG1OmHBKVac7FBerzHX60iSvdJWdQcAGoO9ImeUhovsx/cyJQb1k\n2wh9xI3bKPWn33KdyDzsDwLWm4H5Ue7QLXdAb43u7z90NLDiwpimPV9oYhKv10IPagQB64rGWM19\nqcof8uoz9MC1qd0LXLURuNtErU9aB8m0/UHA+hHkexkS3M99KrV1vYKGrYpySmTWoT/d68RCratn\nGJye4g4C1glNJzdk1RfLfOV6itcsmJ65GKHRgpOihglLrc+zao0wST4tfccMAtYD5WsFImI1RVyC\nOeRTLYZSvQZG3vQmV0yv0FeqR848vH5O32nhP6MbsD7KnqlcJP+VFjFnbOxI9MNl0My33J0GwX6C\nTMvgrtUuM9MjosEDi8VrJttnuZsY/eYoA93FY1oa0Mf8aveEHf8bq7h+AgNBwAJxiRx2TqlHw7Zl\nuDVmc1leL7Zq/firN7FpzAsELAqX3KohG8IG9h/67Sh8l2mvp1SdpV86CkW/Te2IXL2DIghYCF7q\nNWo0Zw4IV5QD+m1cPo8fGHYzOEA/2kVX8y00goD5065eQxyvG/07HtPSUbTv85WdLZewWzcMtTe2\nefglzwsC5oyLem1rytte+ECK6lU8/ah/f6BhEegrAd3XZnbufTNCGr0nmt6Oer2C+zUEjXq1HPDF\n3+Vnix0ZuLwkV/0MrmLwLiBgbmh6u/xtNvXq3M1SXfsOaCKHXwe0b0mUtx4OGeIbO+BXmzHsOAgh\n+tDoWpl7jtzJZ+k52ZR7B1oM0vF5VYUBWSuCG6sJmLCYHJcK1Rg5NPRJZYKyIY8Z+vP1mEZZbDf1\nktHIVfEYNG9tlhKwq6V+7SYJrdRwQG3vqlp1kO9GEjK3rWfn/J4AABUrSURBVAOy0cbdGSE2KJwV\n8bDQGDCz1xqY+7qxMkfLfPrt4JZVh+gVBdsYxMglE5HpoHFrbCeaWSapYZkLmYWlPLAizxTelh1X\n7Vkbynrbt4Vpcpf7062rzzumvD64ludos8kOGxMdU+ez2bmNeY22JxsJ2KtNBw0Eh18oZoreaFvG\n63Zdpwz0qc6Lr/vTWb2qKhWyCjWnt7OMD7TMhYSyi4D5Dpd9sjbcp715FhvomUXcn5TZJl3MVWPM\nSidMcxZswlICdp2z3BbGX7+y0cf36hC00bRB2dTaPBSGm5OnZbpb7BGsXkpx0pfWmHlYPP6rBLRw\nOpYSsOOtF52f9LHL4jw0QhJgar4yD71swGxOyrlarZYIp7jEzdChfdgrC7EdTcfYSr2qWnsmUoY2\nCa5EpzzUnntzLpWnyMfnWS7CvDuDgFXQMhYoO1ieDhC0BpPk6uBHS/AwyFY1atSuWF8l5NFC0ICA\naWnvz+3R/2zU7gSAzrQED4NK/mHQiRYTQpZWBQFT0aheLmvXPfFtML7XEBpFqD2NVkYQledXr9Xd\njrn+V6lYvk5Yqi68CQiYA4uply8LX1p+QvPmx/KV+SIcrxHFr4MhJwhYmaI+mc/VH+NCVQZ/t6Ig\ngjijciy5SieCREXvBRbLaVz/BgMIWIGWtEOv8vvjsnpP5DAt+d0vg1zdTlEmg8gHeOVe3b5NcpMX\nAAFrwmXPciq8ZovMOgcSd/PdTbrFr3JPx38e/6f45WhMfSCrbWT2RRM8bJnJ9lQ4fV3t7teMyr0Y\nS97/2h3Kr8cr38pxPeC4qNT1SM0+buQtFASsicWChy5bf9qvqz3yAwIJDe+VWrnSnPLqhGkqOr+9\nRSmVrZIPmOWJJAQB+6Qlt1CjBE2Nq0QpPO2FHG2/TaMsQTPz3ZaW0XDISGoQKv25ysL1behwf3Da\n9CBg73QYBXoGD5VHmpMtz2NqL+o1MlPkGtVBxraiXbG+DmjR0S9eC4zOC9sKkjgsNE6R+s+wGoMY\nxXzC30J3Va8718bNmYrKNXaYGvN+ZPmAr1McbQmz7AAC9oLG8hqzkFPlbjROPGs9oVPtXJLsmate\nWWw1xdyV5Ij0LS/j+nlLVqSybeAIAvZOy/qWTLYhpjFxo0q9rtJV204B30EHZucrPCgfIxxs4KvX\nZOv+s4OAeZLNOl0SAl3UK0i6TtCwbkTc56rH96pGhlYJZ7WYU204HbttAQG705heWCzcfK6hLpel\nL/lbTRWh0nWChrVTvIdD5mcu4vR6aUUNq11hlad0Xx0q1ax3LhAwN5SG3sdYO6iXshkH/RMacHRl\nDBp2rpZp8kSw9v6QRv8/NLpfSWy3PeW93Q0d0pkjMqEnos+1d0s3r2rDl8ulzGIXGnB+3uibEh6I\nAAFbk8a8jBb1Yh46NRp9mmWioNG545KUGD09rRJO0EAI0Yd2p8exJXGrdPnVi5WwPmS4z7ZUjq/D\nQq8IoQoCAdOSP36oSX+SjxGyLYqFd0vWAIGWEVx/wPWwUTIm2PDtMP26V9AVoV5xIGD/nwXsTNN+\nZc+vKny44wWH6/3Xa9jh8WIar5bryxE07PB+H0dtUiLoYQ1sERrTDqNzOp6nKI809PBZVmgG4rjQ\ndR3xDfe8w5P6uhDh88NjTmYrYXhgdi4QsBWIU6/ivhbh26/jDaegRhE4JmvYBv3Gwfq1efoPhc+P\nhytWZYGa+yBcO9auBwFzYOx8v129hLMOD/WyjQLXU/ColFRJjmbRtMqdfY2/vabefX0VyteFVOXQ\nFwPpxcOqjgEZBExFzgFUKSEDfS8v/+lcsdCPpDkfWSqUkwPDzb/+9yttPQ7BCXs9+FAY6lPJhGOK\nYJxeIGCzouwDZt+rcdHLPfSHK6akVryV62GH9eaHPi+h8YKGfcUS9Wt+LWDDjpCF2A9lZpcGL/Wq\njahoOnloPr3jPVyY2rukPP4vJss8gmJI8Our6KtDvXzBA5uMqpUnzbRaf5ZGupRta4HwoAabH3aU\nHp9Xhl4HOiduaND0oOR3NRt4YDNhyJv4+vZLvcwFdh7UpvADxmLww5Q+1vXI4Q9CvsyvdS+h5Wfw\nwPfqEKcI8MC6YvYequQhInJYTL7q2Tk1KxY4akd9/sVx8UL0q0GhBtDyEAU7KV7mM2vD1ow4fw4O\nBKw/tQNr7ejgGzl0yUWMgMUwJbag3/X2KpUs50gt20lVAsu1nKooSPu8E77YSMDydDClhhka7O57\nCbUb1EuWnOHPJRuOFmuQsVoHy8tluaIc0+UOpfn20C3+XRumaVXxGGhkFwG7WnCGmY4wwzX3f1/1\nKvpe+nnl+bcmB6R42LW0DI8yiAiLtblK1xUjpZz8/riN8oa5jtfDLWpYbY0R+R2+BW7CLgIWjWGI\n+VooNofabZHDr9LkUyJCKLcbsqoyDcdwn23e1fWY14SIZyG1YbrrwY1rorU67Q4Gb2B3AXOx15b1\nGJcwi1COzfdqdLxa7mqV21EbjF0DF4esRZMMA/2rVmkkzZFitHDUFOr5QBez2Dh2F7DZZz2G8P15\nYoR6OfZ8fYRQvw7f3qoM+BqtQclclrvcu54yI+OnT5pF3z4y9tqYxSw2jt0FbGry+F5BCTL6Icl8\nOpyMUrLOaDTsPPII3tccUexW7CJgQtKEY/k9bbGoKLKw6U8pdrPQfritCEVbrFz17w/9ipThlAiU\n1lKVteGSk3JjIr1Pzi4CdqxlK0p/SH9ii3pF39idNSxJA67rVXolG+VhKK3llLGqJKNrLfIBz1Yp\njwQ9GwlYI0nG0BZFMajXKMfrhqxh2ypcN6731qBkr+UEUZtRVczsEE68/vc1J0U4HlxAwFRoekWH\nYbSzerW0BLLhJSS1614GN66R2p7ostZFXxgCAjYNxT7ZJ3JokC6hYVWjDG5WFcKajdcTuZWmVDL9\nKWYM1tI5+RBcQMA8CRpkNT3KN2vDlov4ValX+qJ8e5G3H5pbql9A1S8mKas2n1JLY2yQPIspQMD8\n8R1JNaX5+l6NClGrScd/g0uNca2qlY9Vabc9IVUhQpZClexPnTHfuWHgBQLmjONIqpw/9sna0FyU\nOfaizwcjkCgQcWe+Ao+aaGFVImKQYFzdqcbFLZQsIQiYFv3Q2T7I6vubTb28irqV2T73Z/nBTAdd\nr13Esg39EWtR7TL2bNjtQxgCAhaCeSyuOst33cvcHt+Bpj2cuCH971iVOBkcssNJdSIKtHmlEAEC\nFkVVbzFEJ4rq9fzWpmrFBX/fflvUMKKIN8beDb2S2Vwr95QK3wKvF0WMsT8IWCxCb2mZu5kjh6+q\nZk5XC+qotnVEhE2JbY3WMR3D4Al5JfsYWltb4IFb1hEE7P9THAHNQ6Rv9NycaqH3yeRaihW1gxop\n0T+IxpG6NoPjKElUi4w5GkZEjkbtYiGYQcC60m7EhsihfFbt8cVvO4C8XWlM9jHUohQzjUTZZCw0\n39LXf4rISYETBKyaUUOnpgNUqZE5mb7P5aNSLgStU56FuwQMz4VPZVPNSVIagpa13HNS4EDAarEt\nz7SjyadYSb0cWVgFi5cWfe1Pt8wsY7WrXMVkn3YiIoERUdCdQcAmwDxOOarXj6pe96r0VSW0O2Gj\nJhwZ6Hzhmuy+4vBdJUs9H67vUlmoB7kVCFg1nUNbLer1dXzV58pmvBb1lTZCv43G9rw0tGdwyMN3\nlY/SP8589RSPBkveeWrlCAJmoU+30fQQWxJ87ef6iy222WslgBWyRsypCvoMDmGU10QUM891XBwy\nbLid/4xuALyjVC/hK8elL2We2+9I5cT50E38maga0Hjtx3H8/Zeqwq9nXdMcvo4UDpPNQNmwsRZy\nvUbb6Zh3C3hgRkJnT5qSDZITql6GcdB2ontLFkO+A45ujdILkbMz5EUvpR823JVpiSsMb/zU4IH9\nD1UToojZk3IFe3b1CiJJM/pguO16F7mK0ycTHLLzgKqvjqmeqXwhEAEClgjl+NKnhyhrIUtwFjrc\nZ42MCeFE+SzPhkZiaO1cF5gKQohNeC016wspBoi83K9ie7yCUZoQSh4nLzm2B/08WPi2WI4cTxNC\nx0JvKhpJqkCc17AARRCwVuQgfpEqQ5dHFmHwqjpe/upKn/7J/LQDytzCq/+kWZeSF7cEDdOfkpNZ\n2jk7CJgDtiVc85LvV2lmX8rAREPJPpjdL70pXo8pKpmcpyNrmMHOcyqcfi6YsPH5QcDcuMb35Tns\neXxV+Tb7jsjdoKdNRFxS4nmWXIghMGjwz1y4FY6R5wcBu9MyFbp1afmYKmxLX7aKXMqpgunnEBxX\nMeXShHCi8OgN/pmZ18Yrw6oy2HYoCFgIvvZqVq+g3A164yy0L3PqkSMQcjjx+bltadmmFkKQ83rM\n0eCq0muCII0+Oy3pIX2ilDAdQQ/67+/v6o29HvD86nfW6+f6Qszo16iuCg1JQMBSU+wtnVUKhUtL\nnkfzJUjnt3oNiF70qkW+NPnE4lnZLnYKELAXkljSv9K+5p4BIpgOw/aJ14Ov6Guv1TAvYasqx9ZN\nkowPcLAGlpagHiIUS59cjxZX+7WEqrwGw9rP6/EJZ2MsayUBDywjRd/raHC/zOkbcTAWdEOzR+I0\nv9dcjBONT/blrOgXvfqshBmobcDwBi8JAvbOQGvTjOZD2jZKZpC3J0H3pCrXTs7XuB5WpWGN0UXD\n8TAvq4UQv0IfR2X0YxRVywC1JSyvBJMOWy7xOmVFmrifvsCrhpkDBvNCIHE4SwnY1ZheDSu5qSk7\nQ/Gw5JcZynTXLhttz3xR81h8ylhtuO9VAPSqgH7AXiHEqkyqzoGIsd6Dcm4e0UjGIIFuG49cnoLQ\n1EbjieuMjcUSrhzLXgL23MaRZGeiJmvjeqTt221J8pQNvG48ynktyiWx51mvp1SlwuuLfR6mrMUF\nvdrNa7GdmVjAaneofOU4yVl5mFEomqS41wOqpFo5OehAldEKEbmApjlQjGwre9NXV1UmLi5AHotN\nzsRrYFUPOK1rUpX31aKmXnfA8U52mBxkm3/ob11Pi03bOwBkJvbAnpx7U24L419f6csMaW5N5PDa\nnoiWjAq5dBg3kw/NT8tstNixGPqLS+q8uSVEWaZmYg/sFSHzsGUUiBhEartN8fg+/dAl9Uvpd24y\nsnxlHk6kW0oi8gbJRdyZpTywICL6hsH30rSkTzdunLRWXXv7Ahgsg8Hwop0wTHEsCJgK3zjDAkZ/\nBriqzqoKiNXepQXuKvhSleazibu/GAiYFhcTNy9pdBiday/wXLDRnHVeOBqTkOHDd9pdYr5M0ci5\nWG0NLJTGaPtYFyFo+eH3hzkj/AuhqYwCm/BlsbU7KJTHnypYlRLsvgoAVSBgddhkwGbr/TGLnO+l\nucvhhqya2hDqbJ1BhZYFWugJAlZN1UzNRbqWHIlkdrvebnzZ0qqCZ6DYwWeZj+4AAmbhul/neDNl\n4avkDB/I5Nq/vmXw1SC7L14bJIY8CKHxhuu6dfDXbyEDCJgdIYshp4kre/IoDWNiqyfuGTWu8jYe\nkAr3wDi27Q4C1sqSRlm7oN2OpnszBNTyvGOy8tUuAsl1CVXYzjUzPK4AQZBGnx3fjqdfA5ejKL6g\nXhGYb5dhk99crhUsAx7YOrhHCK8aFiQeJLl4UZsvLh9ZXOX9UbXWO9bNGuuEKadofRqzEggYFIiQ\nsdokl7S7v9emuFd9rtubP5CYuW05QcB2xCUvy6YohtOLeYm1zdiHrwddZQAuo6rjzvRGHeqvYclV\nc2oQsKXwjSO9nnX+/TrunAfI3+pRBkVri90KQcOOLndPGbGsPctQl34N2AvUKxQEbF8ah7DX3W+O\nG+A0bWN0uPI6KRGG7DNZI+4eFh9i/xylnoFE7DMashBXo2qO6Ztq+HehsSj2hDkia1icU/Iv6+ub\n+/hhBLc7gAe2IIYVjjyCoffhmN564R5OVJYmL4xNvTOstnDUzgYClh1zT6s6S5k2HUptTnbjAUsi\nZ20c3/fWxQD6ZNXLT7Z2Gfjo/ipqoT2OzdgEBGxNbMsbzxSMbksFVXW1vPRhGWpnNhqT+MrBUSYQ\nOmaWGpI+DJyy3dn1BC8QsGVp6ZznuXFKZs7I39O10iMIm35aczvAcR9YhyHeJu3n34Ya209H82wg\nYBPQEq9viZNcT3kOYbUFtpdw0NV1yBp21N/Gng7KkEf8nLEdiqtu3BYJ7SBg6+MSJ3mea9h/aq79\nqJy5o3MyQcs/X8ySmKP0O7+Oh/4gYNPQrkC+8cCevZelhVr6ZG0U8U3M0eOSZNhN3TFsMwjYHLjs\nXLkt0U/RbQxN3WdEkIdpjc28Zm203L2gtI7aw44Z3nx47GSrQSBg0+AY9smQNC9jbhspHlf0NiOv\nd1bVWHvKRNMpyAYCNhO25HihtN8f5lmzL15OwG5DoSYzvkokut3AoO0TJ8mdsMxtmwUEbDJODTu8\nd8P8uCViPQ9wxF04NxwRlLFlg4yFYosMGyrqnK6iJGGTJgUBm49ryq97H3gtMCgu59v4DdXrpGpn\n96jQcbuH3bKTJI95oF6OIGCzEuGKCXWFlt9OnuGpP7UDdLfQsWP5LptAhivH8AasBwI2Mc+RaMO+\nwaBwNL857Ie7n+3yUByf7ygZ27l7RoOArYB7MnR+GBRutG9Xz3Yngx7xLY4aUUWSrKgdQMCWQkiG\nXqAX7SPPZty3q/enT/uLiUtehUMoCNiyVAWIUnW5r6amamRaZnTHx062kt8cEEDAdkHupan2/zKg\nuKDcm9zthUnCtzxxsLGdgO2cribAPUmLi8V+leAePTO0AcDMRgKWyskYq6MDa2cCoaeDxc7yksk9\nq4Yi2z2bqzmmkjTowIzWfhtAMdqtmNFie7KRB/YE44DpwGgBTtYUsPUyyGFtsFgAA2sKGP0f5gKL\nBTDwn9ENAAAAsLBdEgcAAKzBmiFEPd1yZL9ekNP5dQmd3zMkVNf/PRHL5EP3uZA8zy6J0Q55s8ky\nRhvEvgLWMx35lrt/s8ghnaFDxyhW1/PC+1QUTbcLSfXs8hjtkHfYg8C+a2B/f39JpjY934aQim4X\nnudZN5LnQvY02s4vLknyrDOzr4Dl4WepGw4H2174Auz57Pa86sxsEULsvMmmqrptJ1nbXriSsUYr\ns+ez2/Oqk7OFgKX9dYZtV2i3vXA9Y41W0LM9n92eV52f3Z/KwCzEX9VJMqy6VXfe8P6/u7jMGDQq\nC3HUsxtb3UCLPRYy2iC4OwAAMCUkcQAAwJQgYAAAMCUIGAAATAkCBgAAU/L/ADTXNdpf43b2AAAA\nAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%matlab\n",
      "% p24.m - pseudospectra of Davies's complex harmonic oscillator\n",
      "%         (For finer, slower plot, change 0:2 to 0:.5.)\n",
      "\n",
      "% Eigenvalues:\n",
      "  N = 70; [D,x] = cheb(N); x = x(2:N);\n",
      "  L = 6; x = L*x; D = D/L;                   % rescale to [-L,L]\n",
      "  A = -D^2; A = A(2:N,2:N) + (1+3i)*diag(x.^2);\n",
      "  clf, plot(eig(A),'.','markersize',14)\n",
      "  axis([0 50 0 40]), drawnow, hold on\n",
      "  \n",
      "% Pseudospectra:\n",
      "  x = 0:2:50; y = 0:2:40; [xx,yy] = meshgrid(x,y); zz = xx+1i*yy;\n",
      "  I = eye(N-1); sigmin = zeros(length(y),length(x));\n",
      "  %h = waitbar(0,'please wait...');                                   %%% removed waitbar commands\n",
      "  for j = 1:length(x),\n",
      "    for i = 1:length(y), sigmin(i,j) = min(svd(zz(i,j)*I-A)); end\n",
      "  end\n",
      "  contour(x,y,sigmin,10.^(-4:.5:-.5)), colormap(1e-6*[1 1 1]);\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAX10lEQVR4nO3d3ZKb\nSrKAUTgx7//KnAtszEZ/SBRVmVVrXUw49tjqaqTmU9II5mVZJgDI5v9aLwAAfiFgAKQkYACkJGAA\npCRgAKQkYACkJGAApCRgAKQkYACkJGAApCRgAKQkYACkJGAApCRgAKQkYACkJGAApCRgAKQkYACk\nJGAApCRgAKQkYACkJGAApCRgAKQkYACkJGAApCRgAKQkYACkJGAApJQ+YPM8b39YtV0PAHX8r/UC\nLtnXa1mWxz8D0KvEE5hQAYws9wT2kSOKAL+JPyFkDdhapv3/vhLhOQgyLEZYRoQ1jLOMM49fag0X\nH+f9P791Q20PXvyr/PaAr/5V2eV9fLQU7/6zHkJc/ppiJApGFuQNQSiVezOmrBPYwbIs2/sFTzNM\nHe3y6oxfKR62lODLOy99wLanoY/nA9LpZm8YlmHulX6+k6d6eqrgvGqv/Ai/xov24EF+l3bx0VLs\nPLP+DgzgCvXqgIBBb4xfPNXf9hQwYDgDjl/91WsSMOiM8auVdAvugIABX0u9s86yeOPXRwIG/ehm\nP5XuG4m84Mhru0jAoBPdHDy81R2Lz/KY/REwgG71HUIBA2K5b5+bZVTquzoFCRj0wPFDHv38ZKW4\nFP0kYEAoxi9vEc4TMEjP+NWNCKfOJ3qWBQyIwvhVSuS1FSRgkJvx66MsK4+wzghrOE/AgBBy7Toj\nrzby2soSMEism/HLwcNSjzlOvSYBAzqWZW8eZJ1BlnGegAEf5B2/7hB5tZHXdgcBg6xG21t9K8v2\ncfDwZwIGvJN0/MrysHGqE2cl5wkYpJRxd/Ooj+8igjG3pIABLyXdLY42fiV9mq4TMMinjx1WH9/F\nz0LVK+lzIWDAc0l3alnGL64TMEimjz1pru/CwcOYBAx44tY9Y67rbhQXrV4pNtpTAgaZ5N3X5GWb\nhyVgwJHx6z7GpoIEDNKwz6ov5jaPcN/LCAQM+A/j131SLDIRAYMcOtj3pfsWYi647KqWZZnneZ7n\nUg9Yk4AB/8TcZX8U+QZd9z1gKcuyrBlrvZCvCRgkEHbfd14H38IV8c/dyDiKCRjwR9LGpBi/Ukg3\nigkYRFdnZ+rcjfseM/74tZeoYf9rvYDfbZt4fUb3W3zAt05ATBmHuSwNSxywaZeu9Q/pXiXwkfGr\n8iOH7U3YhTWU+BDi43OZ7jeQ8J59VmVhDx7yVO4J7JCrw0D25u8AG+PXTQ9YUIW1ZXz3nztg+99+\nvXp2w74i4b0ODh7eJ8WyUyxyc1hqip5lPYR42LgptjUQVuR57spD9f2LlawT2P4kmcNZiIne8sAr\nfYxfbo4cwbq37HKLZQ3Y9BCqLp8eoILOxq/dxPXvDO0uG5b1ECJ0zPhV/2ELir/CbggYMLTIvbm2\ntnmalm0ay/LZ5K8IGMQSeX96nvEriP1a+muYgMGIou1nz4h/3Y0gj7ae1jZNT/5tZw0TMAgkY1ce\n9fFdtFViG86PD9DZU5P4LETgNxn3YkONX9cf6vGyRO8v+JCUgEEUGbvyKH5p7hCqXpvHz8t2RsBg\nLEmvfFhcwOTcMSRleTp+I2AQQqJdf2WDjHQBlxSfkzhgIOnGr/i79SIrjP9txiRgQAGJdsHRlhpt\nPYk4hAjt9XHtqOLij3QXHy3mmYGJPigmYDCEdAcP47ter2gbbR/UFBkTMGgs4I7sK4muGhVnU8dZ\nyZT5XlQCBv0Ltbvsw8WbTAZ5OmIewzxPwKClOPuy34w5fmV/1qb86VoJGHSug73tdXE2QtuV5D1a\n+JSAQTNx9qq/STR+FXTbdz1N/737SenH72HkOhAw6Nl9MQiemb04S321ku2Mv3ku3LDORq4DAYM2\n4uxVQym+WUJ98OvclyjQsL67tREw6Jbxq6wUH/zq8lDhKwIGDYwZgI+Cj191/LbeQUauAwGDPhm/\nEh08vBitabBubQQMassSgMq63yzr7SUvfo+itSdg0CHjV8zx63zDXl2KMMXGr0bAoKoqp7HlaMxe\nzN7c8Wj7hr25YG66Z7AJAYN6MqZlb8xPLhe3NmxSqcsEDLqSLgbBzzy8aXvmeo7CckdmqCT7wcMU\naUyxSEoRMKCZ4L0JvjwEDGowflV4TL0ZjYABPCGH8QkY3M74VeExW5278eZUeO4mYHAvb+QfBa/X\nD1+91ZcenIBBeunGr+C+Gr+WZdk+19WNLN+OgMGNBtz7f9TN+DXP8zT9+bp9NGz+K8uLNvEHmQ+3\nDxjzbgKQa/yKv3M8fa3Cf39Y/3qRa/XWl/rqwIkDNv03Xdumz/gaokteihW028j/xq+MDevjHX/i\ngKXe7hCf8esHkRuWeth6KnHApnO/aTz8nT6eNliF3Vc+uuNXRMVPnf/mb25fN+5T8NXeL+Pv8HIH\n7OMtCSbFohGf/XoUf547+YC7UzaW/T9pfo35K+/XD385Rc+yBizRG09IJ/7Bwztmrx9299so1ipd\ngx9hyhqw/UmrzkIkmuzj11CubMn6d/YavFgHWQM2PTxzgz+RUErwaansAxY5a6PISp56PI5nR7eX\nOGAQU+rxK/5gF6peZcnVtwQMuFG0SGyaL0yurhMwKMn4dasiK2x7rqDf1hckYMA0ZbhnSql6tU2X\nbhUkYFBM6vEruNT1kq6bCBikkevgYbTW1l+Po4V3EzAoI9r++rz4K7++wsrfo5GrDgGDHOJnZi/U\nye41N5101SRgUECuuuz1ffCwWk4cLWxCwCCBm5IQpzSvXFlhne/OyNWQgMFV8TNQU5CDh0//7f6j\nw9fXKF3NCRhEN+z49bOP9br44Osfet16iQgYXHJ3BnJlJsj4dRMjVzQCBiMKmIeY9lcstMWiETD4\nnfFrL874tb9f4H//4/bnzwvY/tXPy+BuAgbDCX7qfJGHevY7sHn356//ecceY5+FgMGPjF+bFEuN\nv8I6Tt7GJUXVBAx+kbdewWMTfHnp9H3XMQEDLol28HBwh2L1vT0FDL5m/LrvAflK3wPWRwIG/Cjg\n/Sq7N9SA9ZGAwXeMX3c8Gm/4LNorAgZfyLvXDr7y4MurT7TOEDAIJMt+3DB3B9H6loDBWXn3s6l7\nk3eznyFaVwgYRNH3nvqV99/1+v92tmVEqxQBg1Py7kPzjl/b19qubZj0KVi5vmJxAgYh5A3kFee/\n6/WvpdtKhq1bCRh8lm6/uYk8fp05eDj9vfDu+hezHE40bNUhYNCeey4f/Lde8zQt8xy9YYat+gQM\nPoi5u6yv/jC31evv//4Rp2Gi1ZaAQWMpxq92wTjWq+liXMkpFgGDd4K80/9W8GV/XN42Y63nHu7/\n4VQ3G4oVmYBBS8FLs2qyyMNxwnme57lSPxwYzELA4KWk1+0NfvDwq1PnK5zON/gdSVITMOhK8Hp9\nq/hXl6ueJA7Y4a2ZqZ+yMo5fzXvzUf3rKB7+S/Dtw1cSB2x6+HC+lyZZxC/NlPb2zd7LjiNxwB5f\nmh1cLY0gUgTmoLODhz+wBxhN4oCt9j9mT6+W5ixYosnYhutu/a6l67rHw63xJQ7Y4SX76rXrNc23\nMgYm+Ph13/mW6x/SPV9BvHl/n6JniQM27TZ3xj0OY4r/Wk1RLyPXbzo7IpU1YOvTsH//5b0YRdx9\npCv+B7/KuiOHkx/zb3R8VkvWgD0+DZ09MXBS5IOHBR/NO9TzOhuz3sgaMLhDxvGroJj1MnKd0fGY\n9YaAQQ3xDx4Wn5auP5p0vTdmtPYEDP6IPyGlsLuK/KUHmUbdKb8nWnsCBrcbZ/ya5z93T56mabuH\n8vePYNf8H6L1ioDBNN05fsUf7Eqv8N/9J883zDkaq6efvhp8m7whYJCgMY8KzkwFv/e/Bw+3vfCH\nRx65W1p1nYDBjYIfPLyt3Ntj/tlLH77KgN1yXfw7CBijyzh+RbZuy78HD/9dK+e/f6fnDW60qkbA\n4C6jjl/TNB1/9dXx7tto1ZCAMbR041e6BfdnqGkyOAFjXOmuuxHwvPnuGbAiEzDIQb3qMGAlImAM\nKt34xX18UjgpAWNE6QJj/CpOtDogYFBY5BtCDl4v0eqMgDGcXDvxXKsNSLQ6JmBQUuTeRF5bWaI1\nCAFjLLl24rlW29yAV6ganIBBMTHvaHzHo8Vh2BqZgDGQRDvxavVKev8twxaTgEEpieq42dacImOG\nLQ4EjFEkCkzlg4drvIp/3VIMW7wiYFBAFx/Vmqdv7qF84zpczIlzBIwhDHvhqPdrW5ZlV4uW34LD\ngw09vYFZCgIGgdQ/83BZlia7L2NWfa9C9XTjp6iagNG/Mcev8wtblmnbWd33rShWTYPcFVrA4Hfd\nfFTrji+rWHUM0qqnBIzOhZ2QDrK30I0fKxi5VU8JGPwobBorLEyu6jDFvidg9Oy+XXnYgenueqX4\nyHNS3hZ8S8DoVpZ6peDTxMXJ1XUCBo1FHr90qyDHA4sTMPqUZfyKWS/dus6AVYGAwRfCHjwstTC/\n4vqZAas+AaNDYTNzULA6Fx/HyPUbV8BqS8DgrJhdvLgqI9dXjFmhJA7Y4T2jt5CsbsrMHedHXH/A\nKw8iXScZs8JKHLDpv+naXlgx3yZDcRfr5cfkDdFKIXHAvKp4NM749fMjGLxeEa10Egdstf4Yv7ny\nv2PW48gyVZQ7d2OavrkOr3TtOc39IMX9Uw4SB+zkT+PgL0quC/jBr/W1P03zNC1n7qE8crq+ugnW\nyA4bJEXPEgds8hJkJ8XBw9KLXBs2vbqZ8oBnNpmrhpI1YOvLdP/zOeDPKpsUBw/vWeTy9JFHGLnc\nW4SsAXt8mXrhUlzYy1s83EN5OfwyuL8fB6MVj7IGDDbBx6+blnd4yLVekbfDt5x+xUcCBs81/6Dx\nt1Lv4g1Y/EDAyC3FuRs8MmBxnYCRWPDMBF9eZYpFcQIGR+kOHobl2hbcSsDIysHDsDo+GZJQBAzK\nG7CChi3qEzBSMn5FIFqdSXH5qD0Bg8K6r6AjhH14/9GFFDETMPKJPH71Wi/DVnZdftJOwGCa+g3P\nRYatpLrM1SMBI5nIpYm8tvMMWxmN+TE7AYNOwnORYSuXMYt1IGBkckdpatbrzM0nq7EHTMENOd8Q\nMKhk3RGt/9tq5+PwYExvTvnzNL0hYKQRf/x6c0+TVuckG7NCMU6VJWDkEL9ecRizmhOqOgSMQd12\nn8nnQ9j+HsrFGbMaetoqT0EdAkYCxWPTZPYq+AUVqwmtikbAiC5dvd78JuwKBwZr0qoUBIzQkv6a\nqkjDjFnVDHLdiv4IGGMJXkRjVgVy1Q0BI67gsXnv/BAmWreSq44JGEH1fd68aN3HoddxCBic9W3/\nDkOYaN3BgDUyASOigOPX+s/X3eVXj+MiuWUZsNgIGOGErdf0d3d5PmN2r9cpFq8IGHzwGL9vM8a3\nHG7lDAEjlmjj15t/K2NliRbfEjACiXOW4Hky9jPHBrlIwOjcTePXwZYxe+E3FIuyBIwooh08/PYL\nTfbIO64lSAUCRghR67WeN7+ef/jyq0wD75rd+IqGBIw+lSviMk1/Gvb4JabB9tQ+NUwoAkZ7MW+Y\nsj8n7vCQg6RLrghOwOhNqRz+9x7Kf+awjtMlV6STPmDb3sqHSJIqOH4Vr8vWsPUqHJ2dZOicQLJL\nHDBvGNm7qS77h0z9AvPzQn8SB+xx8Or48E6vSlWns9moCAMW3UscsEdPP0zqx7h76rXyUudnrz4O\nEVw/AXt/wToCKnSu4Oj1crsWzjv/ub0USeskYPZi6ajXFU5Z4r1BPmDeScC2Ow1O3T1DsBItnhr5\nql3pA7Y9T4M8YX0wfp0kWuyN3Kqn0gcM+uNwAlp1hoBRm/HrFd0alk/p/UbAyKenejlIOCC5KkXA\nqKqn9lxh2BqEI4G3EjDqGfngoU8Z922Q09ajETCSWT8ykWK/4PBgZ958ttfz24SAUUm16mw7mfq7\nFGNWB1QqEQGjhrL1ejOE7Xc+81yjYcasjFSqDwLG7bIc8TtPtLLwq6m+CRj36qZeohWc8/0GJGCk\n9Ooo4nYP5anE78BEKyatYiVg3KjJ+FWwW/aJzWkVbwgYd7m7XmXPpzdsNadVfEvAuEWWX30ZtprQ\nKooQMMqrVq/zQ5iz0Rpy6T9uImCkd+be5/aY1cgV1QgYhVU+eGjn2Jzrj9CKgBGFX0elYMAiDgGj\npJ/Hr/0/XHeRdotBGLAIS8Bo75C99c8GslYUiywEjGKKX7F3e9jJbvRmPgZHRgJGGUUOHj61H8js\nW0sxZtEBAYOe+QAcHRMwrrp7NjJ7neTyFoxGwLjk+u+91qtpPN3TStcrWgWTgHFFkbM21kitDzbt\nriWf5WqKd9MqeEXA+NF9gRl58PIxYThPwPjFDfVapmmepmmeB9plyxVcIWB8rfTnvaZxPrLs5HUo\nSMD4zh1HDnvdjRuw4FYCxhecWPGRS1pANQLGWer1lKOC0IqAcYp6bRQLghAwPlMvBwYhIAHjg2Hr\nJVoQnIDxzoD1ch8yyCJ9wLY9rP1OcePUy7AFGSUO2H6nc7ghvX0QZ3jTA6klDthh8KKsXt8HGLag\nG4kDdpKTnn/QWb1ECz7KOAz0HzA7rB+s9+hKvelEC75y+DFJ0bP+A8YgjNowmk4C9veuvn/+3HYx\nVGPMgpGlD9i227L/Kuv9UcT9tFNzwxuzgE36gFFf5WPjxizgKQEjItECPhIwXnp1FHF3D+WSRAv4\nioDxi61hRUKzpku0gK8IGD+SLqAtAeODmz7RLF3ARQLGO3dc6V+6gCIEjM/2V/o//JevSBdQkIDx\nhZ9LJl1AcQI2kIIVeSzZq0eWLuAmAjaKm25dvX+QwyNLF3ArARvC/kzCUr/QenR4ZOkCbiVg/Xt1\nHnyFkgHc5/9aL4D2lmWRHCAdExiO+AEpmcD6t7/b56P1AKN6AemYwMZl8AJSE7Ah/B3CthMR77rC\nIUA1AjageZ4NXkB6AjaE3a/A/s1hAKk5iWMIyzJN07LVy/QFdEDAhiJcQD8cQhyFqQvojAkMgJQE\nDICUBAyAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyAlAQstHne38oLgH+6uhr9\nvNvZd3DH4e27mWfXkgc46ipgUxfdekrDAA56O4Q4z/PsoBvAAPqcwOZ53kaxQ8+Sjmg5Vw2kkfGt\nf1cBexqnpMWaRAuo6LCrTNGzfg4hptjcAJTSzwS2LMvWsLxTFwAn9ROwSbcARtLPIUQAhiJgAKQk\nYACkJGAApCRgAKQkYACkJGAApCRgAKQkYACkJGDluSgjQAVdXUqquf09lCeXkwe4kwkMgJQEDICU\nBAyAlPwOrKT1l17z7LdfALczgZWnXgAVCBgAKQkYACkJGAApCRgAKQkYACkJGAApCRgAKQkYACkJ\nGAApCRgAKQkYACkJGAApuRr9v9soT67DC5DH6BPYvl4AJDJ6wABIavSAOWYIkNToAdsTM4BEnMSh\nWwApmcAASEnAaphjnOwYYRkR1jBZRrA1TDGWEWENU5hlpNDVIcTtiV8cFgToXT8Bm+d569b+zwB0\nySFEAFLqZ1J5OoE5mgzwm/h16OcQ4lPxnwAAfuMQIgAp9XMIcXIWIsBIugoYAONwCBGAlHo+iSPI\nEcX9mZCtVvJ4Wmb9lRxOE93+e7WVHL73Vpvi1TJqruTxe2+yNSJsim0lbTfFm2XUXMnhiwbZhb7R\nbcBCfa654VffvyJbbZOnH2Zosk32O8qGL49tGesf2m6KdT/VamtE2BQRfkamZz8mDV8YU7Bd6CsO\nIdYwz3OrT6Qty9L8lfd0DfW3SfPtsLIpNhE2RZBd89NlNNlvNNxZ/aDbCSyUw9tMpnbbZJs5an7R\nV8tY/9xkUzTfAqvDMvyk7DXZGrkuAWECu50fxUetjhRF2DMeltFqPRFG88MymvzO6fC/TTxdRtvj\nh1kI2L2yvJGpqeE2CfLz+fQX9dUEeU0eltFkVctfU9PXxuMyRn5hfKX9G9L7BDmFJsIywp6F2OpE\nklZnWEVbRsyzEEd7ccZZRpxn5KSeAwZAxxxCBCAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyAlAQMg\nJQEDICUBAyAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyAl\nAQMgJQEDICUBAyAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyAlAQMgJQEDICUBAyCl/wcMqt45N4Uo\nOgAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    }
   ],
   "metadata": {}
  }
 ]
}