{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Practical quantum mechanics using computers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The standard approach to learning quantum mechanics (QM) involves a fair amount of theory, with some relatively simple pen-and-paper calculations once the theory has been introduced.  But the vast majority of practical QM calculations are done numerically, using a computer.  This series of ipython notebooks is designed to show you how this works, and to give you the ability to do these calculations for yourself.  They should deepen your understanding of how QM works.\n",
    "\n",
    "In this notebook, we will be setting up some basic python functions to allow you to calculate the energies of arbitrary states in an infinite square well; these functions will then be used when we change the basis set, and change the potential.\n",
    "\n",
    "The rough flow of all these QM calculations will be:\n",
    "\n",
    "* Choose and define a basis set\n",
    "* Define the potential in the well (either as a function or by values on a grid)\n",
    "* Create the Hamiltonian matrix\n",
    "* Diagonalise (using libraries) to give eigenvalues and eigenvectors\n",
    "* Display the results\n",
    "\n",
    "In general, if we have a basis set $\\vert\\phi_i\\rangle, i=1\\rightarrow N$, we write the Hamiltonian matrix element between two basis functions as:\n",
    "\n",
    "$$\n",
    "H_{ij} = \\langle \\phi_i\\vert\\hat{H}\\vert\\phi_j\\rangle\n",
    "$$\n",
    "\n",
    "Typically these matrix elements are found by integration on a grid in real space, though there are many other approaches (and some operators and functions where this is not possible).  Remember that the indices $i$ and $j$ in $H_{ij}$ are *arbitrary*: that is, they can each take on any of the values from $1$ to $N$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The eigenbasis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simplest (but least illustrative !) basis set to use is that of the eigenvectors for the system.  For the square well, these are very simple to deduce (they are identical to the solutions for waves on a classical string fixed at both ends).  If we define the potential to be $V(x) = 0, 0 < x < a; V(x) = \\infty, x>a; V(x) = \\infty, x<0$, then the eigenvectors are:\n",
    "\n",
    "$$\n",
    "\\phi_n(x) = A\\sin\\left(\\frac{n\\pi x}{a}\\right), n=1, 2, 3, \\ldots\n",
    "$$\n",
    "\n",
    "where $A$ is an appropriate normalisation constant (it's easy to show that the constant is $\\sqrt{2/a}$ but we want a general approach to normalisation so we'll add a numerical normalisation routine).  We can create a function to generate a set of eigenvectors in an array in python, and another function to normalise them; we will do this in just a moment.  \n",
    "\n",
    "However, we must first note an important point: the basis set we have defined here is technically infinite: that is, when we say that $i=1 \\rightarrow N$, we have $N=\\infty$.  Of course, even with powerful computers, we cannot use an infinite basis set, so we will have to *truncate* the set by choosing a maximum, finite value for $N$ (and note that this will mean that our calculations will be approximations, though we can improve the approximation by increasing the value we have chosen for $N$).\n",
    "\n",
    "In the next cell we will write a function `square_well_eigenfunctions` which will create the general eigenstate $\\phi_{n}(x)$. For completeness (and as a check) we'll then plot some of these.  We also create a function to return the second derivatives (for the kinetic energy), `second_derivative_square_well_eigenfunctions`, and an integration routine that we can use for normalisation, `integrate_functions`.  (In this case, we can write the second derivatives rather easily; we will see a numerical approach to differentiation in another notebook.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0, 11)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD7CAYAAACVMATUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VdW9//H3SkhCCElIyEQmMpCBIUwO1EJrxDpQtahV\nKwgoentb296rXrW1o2hr1au31drqr1YL4oBWrdR5gJJavDhQmTNPZJ5zQuaTnLN+f6wcAlyBDCc5\nQ76v51lPEjjkLPbZ+7PXXnuttZXWGiGEEN7Bx9UVEEII4TwS6kII4UUk1IUQwotIqAshhBeRUBdC\nCC8ioS6EEF5k0li/gVJKxkwKIcQIaK3VcP/NuLTUtdZStObuu+92eR3cpci2kG0h2+LUZaSk+0UI\nIbyIhLoQQngRCfVxlJ2d7eoquA3ZFoNkWwySbTF6ajR9N0N6A6X0WL+HEEJ4G6UU2l1vlAohhBgf\nEupCCOFFJNSFEMKLSKgLIYQXkVAXQggvIqEuhBBeREJdCCG8iIS6EEJ4EQl1IYTwIqcNdaXU00qp\neqXU/mP+LEwp9b5SqkAp9Z5SKnRsqymEEGIohtJS3whcdMKf3QVs01pnAH8HfuzsigkhhBi+04a6\n1non0HrCH68Enhn4/hngcifXSwghxAiMtE89SmtdD6C1rgOinFclIYQQI+Wsx9l5xDKM3TYbhzo7\nOdzbS53VSm1vL0dstqN/76sUkX5+zPD3J8bfn4wpU0iePBmlhr1QmhDDorWmrKeHgq4us29arTT2\n9WE7ZoXTEF9fZgQEEOPvz8yAAOYGBRHo6+vCWgt3NNJQr1dKRWut65VSMUDDqV68YcOGo99nZ2eP\n25rJtb29vNfSwnaLhc/b2ynt6SE9MJCUwMCjwZ3q748jsvu1psFqpbCrixqrlbyuLtr7+1kwdSpL\nQkK4KDycZaGhBPjIoCExOr12O/+0WHi/tZVPjhxhX0cHwZMmMXvKFGIH9s2EgAAmDTQoNNDW38+h\nzk62tbZS2t1NUXc3KZMnsyg4mPOnTeOi8HBmBAS49j8mRiwnJ4ecnJxR/54hraeulEoC3tBaZw38\n/CDQorV+UCn1IyBMa33XSf7tuK6nXtLdzbN1dbzW1ERVby9fCwvja2FhnBUczOygoGEHcqPVyr6O\nDj46coT3Wlo42NnJudOmsSoqissjIpgiLSUxRF02G1ubmtjS0MA/LBbmBQVxUXg4S0NCWDh1KhH+\n/sP6fb12O3mdnXzW3s621lY+aG0lISCAKyMjWRMdTWpg4Bj9T8R4GOl66qcNdaXUC0A2MB2oB+4G\ntgIvAwnAYeAarbXlJP9+zEO922bjhYYGNtbWUtTdzaqoKK6JiuLs4GAmOblV3dLXxzstLTxfX8+u\nI0e4PCKC78bGsiQkxKnvI7zHx21t/LG2lq1NTZwTEsKa6GguDg8n3M/Pqe/Tb7fzaXs7LzU0sKWh\ngfTAQNbPmMHqqCjppvFAYxbqozWWod5gtfJ4dTVP1NRwVnAw342N5aLwcPzGqXuktreX5+vr+UNN\nDbH+/tyekMDKiAh8pQ9+wrNpzdamJn5TWUmN1cr3Y2O5Ljp63LpH+ux23m1p4Y81NXzW3s7NsbF8\nLy6OqGFeDQjXmVCh3mC1cn9FBc/U1XFNZCS3xseTGRTk1PcYjn67na1NTfxPVRXNfX3ck5TEt6Ki\n8JFwn3DsWvNiQwMbysuJ8PPj9oQELnfxiT6/s5NHqqr4S2Mj18fE8OPERAl3DzAhQr2tv5+HKip4\noqaG66Kj+UliIjFudmNoe2srPy0tpctu51fJyVw2fbqMnpkAtNa80dzMz8rKCPL15b7kZJaHhbm6\nWsep7e3l1xUVvFBfz82xsdyZmEjoJGcNgBPO5tWhbtOajbW1/KysjBXTp7MhKYmZkyc7qYbOp7Xm\nzeZmflxaSrS/P4/OmsW8qVNdXS0xRg52dPCfxcXmCjIlhUvd/ERe3t3NPYcP805zM/elpLA+Jkau\nKt2Q14b6x21t/KCoCH8fHx5LS+OM4GAn1m5s9dvt/L+aGu49fJhVUVHck5TENCffHBOuY+nr4xfl\n5bzY0MDdSUl8Z8YMp9+YH0v/am/nP4qK6NOax2bN4kuhsoSTO/G6UG/v7+fHpaX8tamJh1JTWR0V\n5datn1Npslr5aVkZbzY381haGldGRrq6SmKUXm1s5D+Lirh0+nTuS04e9nBEd6G15vn6en5YWso3\nIyP5dXIywdIl4xa8KtTfam7m5sJCLggL4+HUVMK8pHW702Lh3woKmBsUxO/T0mSiiAeq7e3l+0VF\n5HZ28lRGBsumTXN1lZyipa+PO0tK2NbayhPp6Xx9+nRXV2nC84pQb+vv59biYj60WHgyI4Pz3exG\nkzP02Gz86vBh/lRby6OzZnFtdLSrqySGaEt9PbcUF/PvM2bws5kzmeyFY7+3t7by7wUFnDttGo/M\nmkWItNpdxuNDfUdrK+vz81kxfToPpaQw1ct3pt1HjrA2P58FQUE8np7u9Ikownma+/r4fmEh+zs7\n2ZyZyZlePtGso7+fO0pKeLelhU2ZmWR7YePKE3hsqPfa7fy0tJQtDQ08lZHBigl02ddts/GTsjJe\nbmhgU2YmXwsPd3WVxAk+aGlhfX4+V0dF8evk5Ak1M/Od5mb+raCAVVFR/DolBX8PugnsDTwy1Au7\nuliVm0tCQABPZ2YyfYK2Vre1tHB9fj5roqP5ZXKyHDxuwGq387OyMl6or+eZ2bO9sitwKJr7+rgp\nP5/K3l62zJlD+pQprq7ShOFRoa615pm6Ou4sLeXepCS+GxvrsSNbnKXRamV9fj71fX28OGeOLMbk\nQsVdXazKyyPG358/Z2QQ6aEjW5xFa80TNTXcXV7OQykpXB8TM+GP1/HgMaHe0d/PzUVFfN7ezktz\n5siknGNorXmsuppfHj7MH9LSuCZKnj0y3l5qaOAHRUX8YuZMfhAXJ+F1jIMdHXwrN5czgoN5PC3N\n6+97uZpHhPr+jg6uOXSIpaGhPJaWJsvWnsS/2tu55tAhLgoP5zepqV45ysLd9Nhs3FZSwgctLbw0\nd65HTXIbT502G/9RVMSuI0f4y5w5ZEmjbMy4dajb7Xb+XFfHXaWl/CY1lbUxMWP6nt6grb+fbxcU\nUNTdzctz5jBL+jLHTHFXF1fn5pIWGMifMjJkPZQh2FxXx+0lJTyQksKN0h0zJtw61Nfl5rK7vZ1X\n5s5ltgtXU/Q0Wmser6nhnvJynkhP55syE9XpXm1s5ObCQu5OSuJ7cm9nWPI6O7nq0CHOCg7m8fR0\nufJ2spGG+rgMs9DAp2ecIYE+TEopvh8Xx1tZWdxRUsKtRUVY7XZXV8srWO12bi0q4o6SEt7OyuL7\n0n8+bLODgvj0jDOwA0s+/5z8zk5XV8lr/MPyhc8cGpJx636RA2Z0Wvv6uCE/n3qrlb/MnUuiG69S\n6e4O9/RwzaFDxPj7sykz02uWoXAVrTVP1dbyk7IyfjdrFqtklvSI2bTml+XlPFlbS+3Spe7b/TKe\nzyj1ZlprHq6s5OHKSv6cmcklE2iilrO81dzMjfn53JmQwO0JCdLYcKK97e1cnZvLBWFhcoN/BBqs\nVq7Ly6PPbmfLnDnETp4soT5R/NNiYVVuLmuio/lVcrJHLffqKn12Oz8vK+P5hga2zJ7tNQtxuZu2\n/n5uys+nrKeHv8ydK/MthuhDi4XVublcHxPDPUlJTPLxce8bpRLqztdgtbI2L49Om40tc+aQIN0x\nJ1XZ08O1ubkE+/qyefZseZTbGDt2vsXjaWlcLfMtTsquNfdXVPBYVRUbMzOPWyZFQn0CsmvNgxUV\nPFpVxZ8zM2W51C/wVnMzN+Xnc2t8PD9MTJQn/Iyj3UeO8K3cXJlvcRL1Aw2znoHulrgTluKWUJ/A\ndlosrM7L45uRkTyQkkKAdMfQa7dzV2kprzY28oJ0t7iMY75FQVcXW+bMYY6MgAPg/ZYWbsjP58aY\nGDYMdLecSEJ9gmvp6+PfCgoo6+lhy+zZZE7ggye/s5NVeXmkTJ7MnzIyZFljF9Na83RtLT8uK+O+\n5GS+PWPGhL1BbbXb+UlpKS81NrI5M5PzTrFQnIS6QGvNkwMP6P5lUhLfmWCTabTW/LGmhp+Xl/Or\n5GT+fQKHhzvK6+xkVW4uKYGB/DE9fcItlJbb2cnavDziAwJ4OiPjtI9AlFAXR+UN7DyR/v48nZFB\n7AR4bF5Nby83FRTQaLXy7OzZMtHNTfUOjEJ6rr6e/5eezjciIlxdpTFn15pHq6q47/Bh7ktJGXJj\nQ0JdHKfPbudXhw/z/2pq+O2sWazy4Ad3n4rWmi0NDdxWXMzNsbH8dOZM/OSegtv7p8XC9fn5ZE+b\nxm9SU5nmpV1kpd3d3FRQgNVuZ/Ps2cMa4imhLr7Qp0eOcFNBATMDAngiPd2rhj5W9PRwc2EhFb29\nPJ2Rwdle/pg5b9Pe388PS0t5vamJ36elcYUXrW3Ub7fzaHU19x8+zA8TE7k9IQHfYTaqJNTFSVnt\n9qNDHzckJXFzXNywdzB3YtOaJ6qr2VBezi3x8fwoMVGeFuXBPrRY+HZBAfOCgvhdWtr/Gdrnafa2\nt/PtwkKCfX15Mj19xCusSqiL08rr7OS7hYW09ffzWFoaX/HAYX4fWiz8Z1ER0yZN4on0dOk79xI9\nNhv3VVTwRHU1dyQkcFtCgscNzW2yWvl5eTl/bWzk105YklhCXQyJ1pqXGxu5o6SEpaGh3J+cTJIH\nTOUu7+7mx2VlfNTWxsOpqVwdGemV9wgmupLubm4vLuZgZycPp6ayMiLC7T9nq93OH2tq+OXhw1wb\nFcU9SUlOWSROQl0MS6fNxn9XVPD76mpWR0fz08REYtzwsre2t5f7Dh9mS0MDP4iL44eJiQTJzESv\n915LC3eWlBDo48OvU1Lc8sHfNq15rr6ee8rLSQ8M5KHUVKc+CUpCXYxIg9XKAxUVPFNXxw0xMdwS\nH+8Wy/pW9PTwaFUVmwbqdVdi4oQb1zzR2bXmLw0N/KK8nPiAAH6UmMiFYWEub7lb7XZebGjggYoK\npvv58evk5DHpynRJqCulbgNuAuzAAWC91tp6wmsk1D1AVU8Pj1RVsbGujgvDw7ktPp6zgoPH9QDS\nWvNZezu/rari/ZYW1sfEcGt8PPFucJIRrtNvt/NCQwP/U1mJTWtuS0hgdVQUgeN8xdZktfJkbS2/\nr65mXlAQdyQkcMEYnmTGPdSVUrHATiBTa21VSr0EvKW13nzC6yTUPciR/n6erq3lsepqpvj4sC4m\nhuuio8d0REJVTw/PNzTwbF0dXXY7/xEXx00zZhAizwoVx9Bas721ld9UVfHxkSNcGRHBupgYloWG\njtlCbb12O283N7O5vp4dra1cHhHBfyUkMH8cHrjtqlDfBSwE2oHXgEe11ttOeJ2EugfSWvNRWxvP\n1NfzamMjGVOmcFFYGBeHh3NWSMiohkTatOazI0d4t6WF91pbKejq4qrISNZFR7M0NNTll9fC/VX3\n9vJCfT3P1NVh6e/n4vBwLgoP52thYaO+SVnV08N7ra2829LC9tZW5gcFsS4mhqsiI8e1oeGq7pf/\nBO4DuoD3tdZrv+A1EuoertduZ2dbG++2tPBuSwtl3d3MnzqVRVOnMn/qVOIDAojx92eGvz+BxwxD\n67bbqbVaqbNaqertZV9HB3s7Otjf0UFyYCAXh4dzcXg4y0JDPW74mnAPWmsKu7t5b2Df/NBiIS4g\ngIUD+2dqYODRfTPczw9HQtq0prGvj1qrlVqrlYKuLvYM7J9dNhsXDDRgLgwPd9kyG65oqU8DXgWu\nBtqAV4CXtdYvnPA6fffddx/9OTs7m+zs7BG9p3APlr4+9nV2sqe9nYOdndQMBHet1UrPMQ/Gnuzj\nwwx/f2L8/Yn192deUBCLgoNZEBTktdPChWv12e3HBfThnp6jDYuW/v6jr/MBIv38iBnYP2cFBrIo\nOJhFU6eSNHmyS9bdz8nJIScn5+jP99xzz7iH+lXARVrrbw/8vBZYorX+wQmvk5a6EEIM00hb6qO5\n5q0AvqSUmqxMJ+j5QN4ofp8QQohRGnGoa60/xXS57AH2AQp40kn1EkIIMQIy+UgIIdyQK7pfhBBC\nuBkJdSGE8CIS6kII4UUk1IUQwotIqAshhBeRUBdCCC8ioS6EEF5EQl0IIbyIhLoQQngRCXUhhPAi\nEupCCOFFJNSFEMKLSKgLIYQXkVAXQggvIqEuhBBeZPweje1pjhyBsjKoqIDKSqirA4vFlLY26O8H\nux1sNggIgKAgmDoVQkNhxgxTYmNh1iyIiwN5sLIQQ2O1QmmpOf5qaqC2Furrob0dOjuho8Mcez4+\n4OsL/v4wbZopYWEQHw8JCZCYCElJ5u8nEHlIhtYmuD/+GD77DA4ehEOHoLUVkpPNjpGQYEI6LMzs\nOCEh4OdndiofH7MTdnSYHc5iMTthTQ1UV0NxsTkJpKVBVhaceaYpCxeaE4EQE1l9Pezebcq//gW5\nuVBVZY65lBTTMJoxA2JiIDjYNJyCgkyY22wm3Ht6zDFmsUBLi/n3lZVw+LA5DpOTYe5cmD8fliyB\ns882x7GbG+lDMiZmqJeXw7Zt8MEH8OGHJtjPOcd82FlZZgeYOdN5resjR6CwEPbtMzvu7t1m550/\nH847D7Kz4StfgcmTnfN+QrirxkbYvh127ICcHGhoGGzonHEGzJtnwtxZreveXigoMA21vXvhk0/M\nMZiQAMuXw9e+Zo7B0FDnvJ8TSaifit1uPsytW02xWMyH6fhAZ84ENextNzpdXebqICfH7OQHD5qd\n7NJL4bLLICpqfOsjxFg5eBBefx3efNM0ZrKzBxszWVnj3zXZ328aWH//u2nY7doFixfD5Zebkpw8\nvvU5CQn1E2ltWsTPPw8vvQQREeYDW7nSfIDu1sfd1ATvvgtvvAHvvWeuGr71LbjyStPtI4QnKSw0\nx92LL5q+8MsvN42Vr37V3INyJ93dJuBfe82cfOLiYPVqWLXK9M+7iIS6Q20tbNpkis0Ga9aYDyg9\nffzqMFpdXfDWW+ag2LYNLr4YbrwRzj/f9CUK4Y7a280+++c/mxud11wD114LX/qS+zWiTsZmM12y\nzz8Pf/0rLFpkjr1vfnPcu0cndqjb7eYy6o9/NH11V19tPoglS8a/W8XZWlvhhRfMgdLUBP/+7/Dt\nb0v3jHAf+/bB738PL79sulVuvBFWrIBJHj64rqfHdBk99ZTph1+zBr77XcjIGJe3n5ih3tEBmzfD\n735nzqLf/75pGQQHj837udqePfD44/DKK+ZS9pZbzM0lIcZbf79pyT72mBl6ePPNcNNNZpSKNyor\nM+H+1FOm+/aWW+DCC8f0CmRihXp9PTz6KDz5pOmju/VWM3rE01vlQ9XSYlruv/sdpKbCnXeaLhpP\nucQVnqujw+x7v/2t6W++9VZzn8rTW+VD1dNj7hM8+qj5/oc/hOuuG5Ox8BMj1EtL4aGHTL/dqlVw\n++1m+NNE1dcHf/mL2Sb9/fDTn5p+TOl3F87W2moaEb//PZx7Ltxxh+krn6i0Nl29Dz5ohkvedht8\n5ztmHL2TjDTUPaNpV1Ji+unOOgvCwyE/H/7wh4kd6GAmQF13nemWefhhc8DNnm1uEvf3u7p2whs0\nNcFPfmJmRh8+DB99ZLr/JnKgg+kVWL7cjFR7/XUzZDo1FR54wNwwdiH3DvWyMli/3tzwTEgwszPv\nu09uEp5IKdP9snOnuVn8zDMm3J97ztzNF2K4WlvhZz8zNwWbm82Nwj//2bNGkY2XxYvNFfPf/25u\nGjvCvbPTJdVxz1CvqTE3Pc88czDM77lHxmufjlJm9MGOHSbcn3jCTO549VVzuSjE6XR0wK9+ZcK7\nttaE+R//aNZQEac2dy5s2QL/+Ie5ek5LMzeSe3vHtRruFeqtrXDXXSaIAgNNN8u993rEOg1uZ/ly\n03L/zW/M1c2SJaYlIcQXsVpN911ampn1uWsXPP20hPlIzJ5t7vu99ZbpnklPN6P0xumq2T1CvafH\n9Ak7LvX27TM/R0a6umaezdEts3u3uan8ne+YYVh797q6ZsJd2O2mdZmZCW+/De+8Y+ZFzJrl6pp5\nvkWLzDj35583VzuLFpltPNaDU1w6+sWxQ/3kJ2bVwvvvhzlzxrQ+E1pfH/zpT+bq56KLzGV2QoKr\nayVc5R//MMNh7XYzguq881xdI++ltbmhetddZuXJhx82IX8Knjf65cMPTZfAo4/Cs8/C3/4mgT7W\n/Pzge98z63LEx5sT6U9/6vK79WKcFRaatViuv96MM//0Uwn0saaUGc9/4ICZ8f71r5vtX1Xl9Lca\nVagrpUKVUi8rpfKUUoeUUktO+4+Ki80iVevWwX/9l1mp8KtfHU01xHCFhJh+9n37zE6VkWFmyslI\nGe/W0mJC/MtfNiU/36yLJJPWxs+kSWapgYIC07BasAB+8Qtzg9pJRvtpPgq8rbWeDSwA8k76Sotl\ncMLCWWdBXp6ZQCQ7lOvEx5vhj6+/bq6WFi+Wm6neqK/PTBzKzDQ3RHNzzUxIWb/fdRwNqz17zDyc\nzEwzv8RuH/WvHnGfulIqBNijtU49zeu0fvxxMyTxsstMP2509IjeU4whrc3So3feaUYfPfSQGQkh\nPJfW5sbc7bebUSy/+Y10cbqrTz4xs1KtVrMEw1e+Mv7LBCilFgBPArmYVvpu4BatdfcJr9M6Oxse\necRcagj31ttrWnUPPmj6/H7+cxlS6okOHTLdmxUVJsxXrHB1jcTpaG2GQv7oR7BkCerll8f9Rukk\nYDHwB631YqALuOuLXrjhq19lw2uvsWHDBnJyckbxlmLMBQSY1vqhQ6afLzPTTGKSZQc8Q1OTmbh3\n3nnmKVr790uge4icf/yDDfn5bFizhg1NTSP+PaNpqUcDu7TWKQM/LwN+pLW+7ITXuf5xdmLk9u0z\nl4UNDabFd+GFrq6R+CKOyUP332/uVW3YYNZJEh5r3Ic0aq3rgUqllGMxiPMxXTHCmyxYYJ6het99\npgV4ySVm1IRwD1qb4cBz55rP6cMPTfeZBPqENarJRwP96k8BfkApsF5r3XbCa6Sl7i2ObQ1ec41p\nDcqsX9f517/MTdCmJvif/zETyoTXcMnkI631Pq31WVrrhVrrK08MdOFl/P3Nzbe8PLNm++zZ8N//\nbZZ5EOOnstLcxL70UjPOfO9eCXRxlAwSF8MXEWEu8f/3f83CTxkZZpy7E8bYilOwWMw084ULzRyD\nwkLzzNqJ8tQhMSQS6mLk0tPN2PYXXjAjZBYvNgtCSXebc/X0mJvUGRmmq2X/fnOPw1ufxStGxbMe\nZyfcl9awdatZSyYiwvS7L13q6lp5tv5+M8vw3nvNCfO++8wNUTEhTIxnlAr3Z7OZrpgNG0wA3Xsv\nnHGGq2vlWWw2MwnlnnsgLg5+/Wt5fNwEJKEu3Etvr1kk7P77TStzwwbzVZyczQYvv2xOhGFhJtTP\nP9+s8CcmHAl14Z56eswa7g88MLjU75e/7OpauZe+PvMghQceMEsy3HsvXHCBhPkEJ6Eu3FtPj1kR\n8sEHITHRrG9x0UUTe5XOjg7YuNE8MCEtzTws5rzzJMwFIKEuPEV/v+kvfvhh00Vz222wZo15Ju1E\nUVVlJnE99RRkZ5u1dpac/lEEYmLxvCcfiYlp0iS47jr4/HP4wx/MWu6JiWat/eJiV9du7GhtpvFf\ndZVZ2ri72zxx6JVXJNCFU0lLXbhecbF5MO+mTea5jTfeaB795Q2t99paeO45ePppc0K7+WZYu9Y8\nJEGIU5DuF+H5enrg1VdN3/vu3aZVe911sGyZWZbAU7S3wxtvmKGdH39sHt+4fr0Zty/95WKIJNSF\nd6mqMiNCXnoJamrMg5KvvNI8z9YdH8PW2Ghm0776KuTkmBPR6tVwxRUwZYqrayc8kIS68F4lJSYs\nt26FgwdNsF98MZx7rpng5IoRNI4+8e3b4d13zYOEzzsPvvlN89hGeVqUGCUJdTExtLTAtm0mSP/5\nT2huNuPezz7bjINfsMDceHVmN0dfn1lDfu9eU3btMg8PmTfPjF65+GLTteLv77z3FBOehLqYmOrq\nYOdOs7b43r0mbDs7ISVlsMTEmPVoIiPNDcqAAFN8fU1g9/aa0tJiulEaG6G6GkpLTamogJkzzQlj\nwQI45xxzEgkKcvX/XngxCXUhHFpaBgO5tNQ8is8R1h0dgyHe329a146QDwszwR8RATNmQGqqOSkk\nJ3vHSBzhUSTUhRDCi8jkIyGEEBLqQgjhTSTUhRDCi0ioCyGEF5FQF0IILyKhLoQQXkRCXQghvIiE\nuhBCeBEJdSGE8CIS6kII4UUk1IUQwotIqAshhBeRUBdCCC8ioS6EEF5k1KGulPJRSn2ulHrdGRUS\nQggxcs5oqd8C5Drh9wghhBilUYW6Uioe+DrwlHOqI4QQYjRG21L/LXAnII82EkIINzBppP9QKXUJ\nUK+13quUygac+Ph279Vh7aD6SDW1HbU0djbS0NlAY1cjlh4Lbb1ttPW00WHtoLu/m66+Lnr7e+m3\n92PTNuzajo/ywVf54qN8CJgUQOCkQAL9AgnyCyIkIITQgFCmTZ5GxJQIIoMiiQqKIjoomriQOKYH\nTkcp+ZiE8GYjDnVgKfANpdTXgUAgWCm1WWu97sQXbtiw4ej32dnZZGdnj+Jt3ZvNbqPcUk5hcyGF\nzYWUtpZSZimjzFJGRVsFVpuVuOA4ZgTPIDoomsgpkURMiSAxNJFpk6cRGhDKVP+pBPoFMsVvCgG+\nAUzymYSvjwlyu7Zjs9uwaRtWm5Xuvm66+7vptHYePSlYeiyUW8r5rOYzGjobqOuoo7q9mu6+buJD\n4kmalkTytGSSw5JJC08jIyKDWeGzmOI3xdWbT4gJKycnh5ycnFH/Hqc8eFopdS5wu9b6G1/wd175\n4GmtNbUdteyp3cP++v0cbDzIgfoDFLUUERUURcb0DNKnp5MalkpyWDLJ05JJCE0gbHKYy1rLXX1d\nVB2poqy1jHJLOaWtpRS1FFHQXEBpaykxU2PIispiXtQ8sqKyWDRjEWnhafj6+LqkvkJMZCN98LSE\n+hA1dTX6U0vkAAAVaUlEQVTxafWnfFL1CZ/VfMbntZ9j0zYWxSxiQfQCE4TRWWRGZHpki7ff3k9p\naymHGg5xsOEg+xv2s6d2D/Wd9WRFZXFW7FksiV/CkrglpISlSDeOEGPMpaF+yjfwwFDXWlPSWsI/\nD/+TnRU72Vm5k7qOOs6MPZMlcUs4O+5szphxBvEh8V4fbpYeC3tq9/BZzWd8Uv0Jn1Z/Sm9/L0sT\nl7IsYRnLEpexeMZi/Hz9XF1VIbyKhPoolVvK2Va6jZzyHHLKcwA4N+lcliYsZVniMuZGzpVuiAGV\nbZV8VPkROyt28s+Kf1JuKefLCV8me2Y256ecz6KYRbKthBglCfVh6rB2sL10O++VvMcHpR/Q1tPG\n+SnnszxpOecln0dqWKrXt8KdpbmrmQ8Pf8iO8h1sK91GQ2cDy5OXc0HKBaxIW0F8SLyrqyiEx5FQ\nH4LC5kLeKHiDt4vf5tPqT1kSt4SLUi/iwtQLyYrOwkfJUjjOUHWkim2l23iv5D3eL3mf2OBYVsxa\nwWXpl3FOwjlM8hnNoCshJgYJ9S9gs9vYVbWLrflbeaPwDTqsHVyadimXpF/C8uTlTPWf6pJ6TSQ2\nu41Pqz/l7aK3ebPoTSrbKlmRtoKVGSu5eNbF8hkIcRIS6gOsNivbS7fz17y/8nrh68RMjeHyjMv5\nRsY3WDxjsXSpuFhlWyVvFr7J1oKt7KrcxblJ53JF5hWszFjJ9CnTXV09IdzGhA713v5e3i95n1fy\nXuGNgjeYHTmbKzOv5IrZV5ASljKm7y1GztJj4e2it3k171W2lW5jSdwSrppzFVfOvpKIKRGurp4Q\nLjXhQr3P1sf2su28dOgl/pb/N+ZFzePqOVdz5ewriQuJc/r7ibHVae3kneJ3eDn3Zd4tfpcvxX+J\nb839FldkXkFYYJirqyfEuJsQoW7Xdj6q+IgtB7fwSu4rpISlcO28a7l6ztUS5F6k09rJW0Vv8dKh\nl9hWuo3spGxWzVvFZemXEeQf5OrqCTEuvDrUDzUc4rn9z/H8gecJCQhhddZqrp13rXStTABtPW1s\nzd/KCwdf4JOqT/hGxjdYM38Ny5OXyyga4dW8LtTrOup44cALPLv/WRo7G1mdtZo189cwP3r+GNRS\neIK6jjpeOvgSzx14jqojVayat4p1C9axIHqB3AAXXscrQr27r5u/FfyNZ/Y9w67KXVyeeTlr568l\nOylbZiiK4xQ0FfDs/md5dv+zhASEsG7+OtbMX8OM4BmurpoQTuGxoa61ZlfVLjbt3cQrua9wZuyZ\nrFuwjisyr5D+U3Fadm3nw8MfsnnfZl7Lf41z4s/h+gXXszJzJZMnTXZ19YQYMY8L9aojVTy771k2\n7duEQnH9gutZu2CtTCkXI9Zp7eS1/NfYtHcTe+r28K253+KGhTdwVuxZ0j0jPI5HhHpPfw9/y/8b\nG/du5NPqT7l6ztWsX7SeJXFL5KATTlXRVsHmfZvZtHcTAZMCWL9wPWvnryV6arSrqybEkLh1qO+u\n3s3GvRt58eCLLJqxiPUL13NF5hUE+gWO6XsLobVmZ8VONu7dyGv5r/HVmV9l/cL1XJJ2iSwXLNya\nW4d68iPJ3LDwBq5fcD0zp80c0/cT4mQ6rB28fOhlNu7dSEFzAWuy1rB+0XrmRc1zddWEOKq4pZhN\nezdx3/n3uW+o2+w2WQFRuJWi5iI27d3EM/ueYUbwDG5ceCPXzrtWZq8Kl+iwdvBK7its3LuRvMY8\nrsu6jkdWPOK+oe4uS+8KcSKb3cYHpR+wce9G3it+jxVpK1i/cD3nJ58vw2jFmDqxa3BZ4jJuXHgj\nl6Rfgr+vv3t3v0ioC0/Q0t3ClgNb2Lh3I/Wd9aybv44bFt5A2vQ0V1dNeJGKtoqjI//8fPy4cdGN\nrJm/hpipMce9TkJdCCc6UH+ATXs38dyB55gVPosbFtzANXOvIXRyqKurJjyQY7jtM/ue4fPaz7lm\nzjWsX7T+lMNtJdSFGAN9tj7eKX6HZ/Y9w/bS7axIW8G6+eu4IPUCWXtGnNJoJ8ZJqAsxxpq7mnnp\n0Ets3reZw22HWTVvFWvnr2VhzEKZZyGOym/K59l9z/LcgecIDQhl3YJ1XJd13bCXsJBQF2IcOdae\neW7/cwT5B7Emaw2rs1bLkN0Jqq6jjhcPvshz+5+jpr1mcLG5mAUj/p0S6kK4gF3b+d/K/+X5/c/z\ncu7LZEZksjprNVfPuZrIoEhXV0+MIUuPhb/m/ZUtB7ewu2Y3KzNWsmb+Gs5LOs8pI6ck1IVwMavN\nyvsl77Pl4BbeKnyLJfFL5OlNXqa9t503Ct/gL4f+wo7yHZyffD6rs1ZzSdolTp8hL6EuhBs58elN\nyxKXcdXsq1iZuZLwwHBXV08MQ3tvO28Wvskrea8c/SyvmXMNKzNXMm3ytDF7Xwl1IdzUsaHwQckH\nnJNwDldkXsHKjJWy/rubaupq4s3CN3kt/zV2lO3gKzO/Mu4nZQl1ITxAh7WDd4reYWvBVt4uepvM\niExWZqzksvTLmBM5R0bRuFBxSzFvFLzB64Wv83nt53wt5WtcnnE5l2VcNqYt8pORUBfCw1htVnLK\nc3i94HXeKHwDX+XLpemXsmLWCrKTsmUV0zFmtVn5qOIj3il+hzcL36S1p5VL0y7l0vRLuTD1Qpdv\nfwl1ITyY1poDDQd4q/At3il+hz11e1iWuIwLUy7kgtQLmBs5V1rxo6S1prilmA9KP+CD0g/YUbaD\njIgMVsxawSVpl3BG7BlutfCghLoQXsTSY2Fb6TY+KDEB1NPfw/Lk5WQnZZOdlE1qWKqE/BBUtlWy\no3wHOeU5/L3s7/TZ+7gg5QIuSLmAC1MvdOthpxLqQnixkpaSo+G0o3wHCsWyxGUsS1zG0oSlZEVn\nTfhlC+zaTl5jHh9VfsTOip3srNhJu7Wd7KRszks6j+ykbGZHzPaYk+G4h7pSKh7YDEQDduBPWuvf\nfcHrJNSFcCKtNSWtJXxUMRBelTupOlLFophFLIlbwpmxZ7J4xmJSw1PdqjvBmbTWlFvK2VO3h901\nu/mk+hN21+wmYkqEOdklLGNp4lKPCvETuSLUY4AYrfVepdRU4F/ASq11/gmvk1AXYoxZeix8Vv0Z\nn1R/wr9q/8We2j20dLcwP3o+WVFZzIuax7yoeWRGZBIVFOVRQdfc1UxBcwEHGw5yoP4ABxsPsq9u\nH4F+gSyesZjFMYtZEr+Es+POJmJKhKur6zQu735RSm0FHtNabz/hzyXUhXCB5q5m9tXv41DDIQ40\nHOBgw0Hym/KxazsZERmkhqWSPC2Z5LBkkqYlER8ST1xwHMEBweNaz66+LqqPVFN1pIrDbYcpay2j\nzFJGSWsJBU0F9Nn7yJiecfTENC9qHguiF3j9Q8RdGupKqSQgB5inte444e8k1IVwI01dTRQ0FVDS\nWnI0QMst5VS3V1N9pJpJPpOInhpN5JRIooKiCA8MJzQglGmTpxESEMIUvykE+gUSOCmQST6T8PXx\nxVf5opTCZrdh0zb67f1093XT3d9NV18XR3qP0NbTRltvGy3dLTR0NtDY1UhDZwPdfd3EhcQRFxzH\nzGkzzYlmWjIpYSlkRGQQHRTtUVcWzjLSUB/1nZWBrpdXgFtODHSHDRs2HP0+Ozub7Ozs0b6tEGKE\nIqZEEJEYwdLEpf/n77TWWHosx4VuS3fL0UAut5TT3W/Curuvm357PzZtw2a3odH4Kl98fXyZ5DOJ\nwEkm+AP9AgkJCCFmagyZEZmEBYYRFRR13EljIob2iXJycsjJyRn17xlVS10pNQl4E3hHa/3oSV4j\nLXUhhBgml3S/KKU2A01a6/86xWsk1IUQYphcMfplKfAhcADQA+UnWut3T3idhLoQQgyTy0e/nPQN\nJNSFEGLYRhrq3jkzQQghJigJdSGE8CIS6kII4UUk1IUQwotIqAshhBeRUBdCCC8ioS6EEF5EQl0I\nIbyIhLoQQngRCXUhhPAiEupCCOFFJNSFEMKLSKgLIYQXkVAXQggvIqEuhBBeZNTPKBViuLSGnh44\ncgTa2szX9nZTOjpM6ewcLF1dg6WnB7q7zdfe3sFitUJfn/na3z9YbDZT7PbBorUpDkqZ4uMDvr6D\nX319YdIkU/z8TPH3NyUgwJTJkwdLYCBMmTJYgoJMmTp1sAQHmxIaCiEhpvj5ue6zEN5HHpIhRkxr\nsFigqWmwNDcPlpYWaG0dLBaLCXGLxQRmSIgJN0fQBQcPBp8jEIOCBkMyMNAUR4g6gjUgwAStI3Qd\nIewIZkfx8RkMcEdxBLzWxwe/44TgODn09Q2eNKzWwZNJT89g6e42patr8GTkOEF1dAyeuNrbzYnM\nUfz9zXaYNs2UsLDBMn26KeHhEBExWCIjzTYR3kuefCScQmsTyLW1UFc3+LW+frA0NEBjownxwEAT\nMBERJnwcX8PDB4sjoByhFRpqgliY7d3VNXjCc5wAW1oGT4rHnigdJ8/GRnOSioyEqChToqMHy4wZ\nEBNjSmysOVkKzyKhLk6rtxeqq6GqarBUVw+WmhoT4EFBJhSODYZjA8MRIhEREs6uorW5AmhoGCzH\nnnjr6gZPyjU15kolNtaUuLjBEh8/WGJizIlCuAcJ9QlOa9OqKy+Hw4cHS0WFKZWV5u9jYiAhwRzE\njoPacYDHxpognzzZ1f8b4Uxam26emhpTjj2RO07ulZXmaiE21uwfiYkwc+bg16Qk8710+YwfCfUJ\noKsLysqgtNSUsrLBUl5u+oiTk81BeOwB6ThIo6NNi02IL9LTY4K+stI0BI5tGJSXm6+hoWYfS0oy\nX1NSTElNNQ2ESTL0wmkk1L2ExQLFxceXkhJTWlrMweQ4kJKTB0tSkum3FmKs2O2mS+fYxoSjkVFS\nYrqAEhJMwM+adXxJSZGuuuGSUPcgXV1QVGRKYeFgKSoyraVjDwbHAZKaarpIpM9TuKueHtOiLyk5\nvjFSVGRa+TExkJYG6enHl5kzpYX/RSTU3YzW5iZVfv7xpaDAtGiSkyEjY3DHduzs0dGmG0UIb9Lf\nbwLf0ZBxfHUcDykp5njIzBwsGRmmu2eiklB3EZvNXILm5kJenvnqCPCAAJg9e3AHdXyVlokQgxxX\nrgUFphzbCAoNNceQo8yZY75OhMaPhPoYs9nMpeShQya4HV8LC83wPsfO5iiZmWa8thBiZOx2MzIn\nL2+w5OaaojXMnWuOuzlzzPdz55ouHm8Jewl1J9Ha9P8dPHh8KSgwrQPHjuTYiTIzZWKHEONJazP5\n6tgGlqM4wn7evMGSleWZgwgk1EegpQUOHDi+HDpkQvrYncIR5BLeQrgvrc3EK0fAO47ngwfN0hNZ\nWceX2bPde06GhPop9PWZlva+fbB/vykHDpgJGY4zuaPMmyfdJkJ4E8fV97GNt/37TXdqcrI57ufP\nN2XBAjMs0x26cCTUBzQ2Dob3vn2mFBaayTeOD27+fPNBzpwpQwSFmKh6e01jz9HQc5SursGAX7DA\nfD9v3vjPpp1woW6zmbB2BPfeveZrV9fgB+H4UObOlenNQoihaWw8vlG4b58J/5kzTZ4sXDj4dcaM\nsWvVuyTUlVIXA49gHrbxtNb6wS94zahDvaNjcCPv3WvKwYPmTrdjAzs2cmKie1w6CSG8R1+fGWJ5\nbANyzx5zpX9syC9caIYtO2PI8riHulLKBygEzgdqgM+Aa7XW+Se8blihXldnNtbevYNfq6pMa/vY\njTd/vlmP25Pk5OSQnZ3t6mq4BdkWg2RbDPKkbeGYYOhoaDoCv6rKDKxwhPyiRSavhjvQYqShPprz\nydlAkdb68EAFXgRWAvmn/FcD7HYzlfjEAO/rG9wQ3/gG/OIXZtigN0zW8aQddqzJthgk22KQJ20L\npQaXM/761wf/3NGz4Mi1jRvN0MuEhOODftEiM0za2UYTlXFA5TE/V2GC/v/o7TXdJY7/5J495j8d\nEWH+YwsXwve+Z76Pj5fuEyGE55o6Fb78ZVMcHCPwHI3Xhx4y3wcEDAa8I+xTUkY3gGNc2r/TpplF\nqRyVv+oq043iiRMChBBiuPz8Bue9rF1r/kxrs8yxI+iffx7uuMOs1LpgwcjfazR96l8CNmitLx74\n+S5An3izVCnlnjOPhBDCzY33jVJfoABzo7QW+BRYpbXOG9EvFEIIMWoj7n7RWtuUUj8A3mdwSKME\nuhBCuNCYTz4SQggxfpw2SV4pdbFSKl8pVaiU+tFJXvM7pVSRUmqvUmqhs97b3ZxuWyilViul9g2U\nnUqpLFfUc6wNZZ8YeN1ZSqk+pdSV41m/8TTE4yNbKbVHKXVQKbVjvOs4XoZwfIQopV4fyIkDSqkb\nXFDNcaGUelopVa+U2n+K1wwvN7XWoy6Yk0MxMBPwA/YCmSe8ZgXw1sD3S4CPnfHe7laGuC2+BIQO\nfH+xN26LoWyHY163HXgTuNLV9XbhPhEKHALiBn6OcHW9Xbgtfgzc79gOQDMwydV1H6PtsQxYCOw/\nyd8POzed1VI/OhFJa90HOCYiHWslsBlAa/0JEKqUGoOh9y532m2htf5Ya9028OPHmDH/3mYo+wTA\nfwCvAA3jWblxNpRtsRp4VWtdDaC1bhrnOo6XoWwLDQQPfB8MNGut+8exjuNGa70TaD3FS4adm84K\n9S+aiHRiUJ34muoveI03GMq2ONa/Ae+MaY1c47TbQSkVC1yutX4C8OYpZ0PZJ9KBcKXUDqXUZ0qp\nteNWu/E1lG3xe2COUqoG2AfcMk51c0fDzk0vmHzvuZRS5wHrMZdgE9EjwLF9qt4c7KczCVgMLAeC\ngF1KqV1a62LXVsslLgL2aK2XK6VSgQ+UUvO11h2urpgncFaoVwOJx/wcP/BnJ74m4TSv8QZD2RYo\npeYDTwIXa61PdfnlqYayHc4EXlRKKUzf6QqlVJ/W+vVxquN4Gcq2qAKatNY9QI9S6kNgAab/2ZsM\nZVusB+4H0FqXKKXKgExg97jU0L0MOzed1f3yGTBLKTVTKeUPXAuceGC+DqyDo7NRLVrreie9vzs5\n7bZQSiUCrwJrtdYlLqjjeDjtdtBapwyUZEy/+ve8MNBhaMfH34BlSilfpdQUzE0xb5z3MZRtcRj4\nGsBA/3E6UDqutRxfipNfpQ47N53SUtcnmYiklPqO+Wv9pNb6baXU15VSxUAn5mzsdYayLYCfA+HA\n4wOt1D6t9Rcuhuaphrgdjvsn417JcTLE4yNfKfUesB+wAU9qrXNdWO0xMcT94lfApmOG+f1Qa93i\noiqPKaXUC0A2MF0pVQHcDfgzityUyUdCCOFF5AmdQgjhRSTUhRDCi0ioCyGEF5FQF0IILyKhLoQQ\nXkRCXQghvIiEuhBCeBEJdSGE8CL/H+RvKJ6jzY0qAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108e68210>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Import libraries and set up in-line plotting.\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as pl\n",
    "import numpy as np\n",
    "\n",
    "# Define the eigenbasis - normalisation needed elsewhere\n",
    "def square_well_eigenfunctions(n,width,norm,x):\n",
    "    \"\"\"The eigenbasis for a square well, running from 0 to a (width), sin(n pi x/a). \n",
    "    N.B. requires a normalisation factor, norm.\"\"\"\n",
    "    wavevector = np.pi*n/width\n",
    "    return norm*np.sin(wavevector*x)\n",
    "\n",
    "# We will also define the second derivative for kinetic energy (KE)\n",
    "def second_derivative_square_well_eigenfunctions(n,width,norm,x):\n",
    "    \"\"\"The second derivative of the eigenbasis for a square well, running from 0 to a, sin(n pi x/a)\"\"\"\n",
    "    wavevector = np.pi*n/width\n",
    "    return -wavevector*wavevector*norm*np.sin(wavevector*x)\n",
    "\n",
    "# Define the x-axis\n",
    "square_well_width = 1.0\n",
    "number_of_x_points = 101\n",
    "x = np.linspace(0.0,square_well_width,number_of_x_points)\n",
    "x_spacing = square_well_width/(number_of_x_points - 1)\n",
    "\n",
    "# Integrate two functions over the width of the well\n",
    "# NB this is a VERY simple integration routine: there are much better ways\n",
    "def integrate_functions(function1,function2,number_of_points,dx):\n",
    "    \"\"\"Integrate the product of two functions over defined x range with spacing dx\"\"\"\n",
    "    integral = 0.0\n",
    "    for i in range(number_of_points):\n",
    "        integral = integral + function1[i]*function2[i]\n",
    "    integral = integral*dx\n",
    "    return integral\n",
    "\n",
    "\n",
    "# Now plot the first few functions; we offset them vertically by adding 3(n-1) \n",
    "# to make it clearer; remember that range(1,5) will give number 1 to 4. We don't\n",
    "# normalise in this plot.\n",
    "for m in range(1,5):\n",
    "    pl.plot(x,square_well_eigenfunctions(m,square_well_width,1.0,x)+3*(m-1))\n",
    "# Set y-axis so that we can see all the functions\n",
    "pl.ylim((0,11))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do these make sense ? You should ask yourself if they match the boundary conditions, and if they fit with your expectations of the physical properties of the system.\n",
    "\n",
    "We should now be confident that our function is correct (it's always worth checking this somehow, whether through plotting the output or a simple test).  So we will create arrays to store the first ten basis states, properly normalised.  I have chosen ten for no particular reason - we could use fewer or more states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# These arrays will each hold an array of functions\n",
    "basis_functions = []\n",
    "second_derivative_basis_functions = []\n",
    "\n",
    "number_of_basis_functions = 10\n",
    "# Loop over first num_basis basis states, normalise and create an array\n",
    "# NB the basis_array will start from 0\n",
    "for n in range(1,number_of_basis_functions+1):\n",
    "    # Calculate A = <phi_n|phi_n>\n",
    "    eigenfunction_n_unnormalised = square_well_eigenfunctions(n,square_well_width,1.0,x)\n",
    "    integral = integrate_functions(eigenfunction_n_unnormalised,eigenfunction_n_unnormalised,number_of_x_points,x_spacing)\n",
    "    # Use 1/sqrt{A} as normalisation constant\n",
    "    normalisation = 1.0/np.sqrt(integral)\n",
    "    basis_functions.append(square_well_eigenfunctions(n,square_well_width,normalisation,x))\n",
    "    second_derivative_basis_functions.append(second_derivative_square_well_eigenfunctions(n,square_well_width,normalisation,x))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how we have done the numerical integration: it is a *very* simple integration scheme (there are far better ones available) which is almost an implementation of the trapezium rule.  How do we calculate the integral with the Hamiltonian ? In this case we can construct the kinetic energy explicitly, and the potential can be specified by its value on the grid.  If we do not have the second derivatives analytically, or do not want to code them, we might also use a numerical approach to differentiate (finite differences).\n",
    "\n",
    "We know that the energies for the eigenvectors are given by $E_{n} = \\hbar^{2} k_{n}^{2}/(2m) = \\hbar^{2} n^{2} \\pi^{2}/(2m a^{2})$ , for a well with width $a$.  We can check this easily in our code:\n",
    "\n",
    "* Calculate $\\langle \\phi_{n} \\vert \\hat{H} \\vert \\phi_{n} \\rangle$ where $\\hat{H} = -\\frac{\\hbar^{2}}{2m} \\frac{d^{2}}{dx^{2}}$ in position representation\n",
    "* For simplicity, we'll set $\\hbar = m = 1$ (atomic units) but we will define the constants in the code\n",
    "* We'll calculate this and check it against the analytical result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy calculated for eigenstate n= 1  is  4.93480220054\n",
      "Energy expected   for eigenstate n= 1  is  4.93480220054\n",
      "Energy calculated for eigenstate n= 2  is  19.7392088022\n",
      "Energy expected   for eigenstate n= 2  is  19.7392088022\n",
      "Energy calculated for eigenstate n= 3  is  44.4132198049\n",
      "Energy expected   for eigenstate n= 3  is  44.4132198049\n",
      "Energy calculated for eigenstate n= 4  is  78.9568352087\n",
      "Energy expected   for eigenstate n= 4  is  78.9568352087\n",
      "Energy calculated for eigenstate n= 5  is  123.370055014\n",
      "Energy expected   for eigenstate n= 5  is  123.370055014\n",
      "Energy calculated for eigenstate n= 6  is  177.65287922\n",
      "Energy expected   for eigenstate n= 6  is  177.65287922\n",
      "Energy calculated for eigenstate n= 7  is  241.805307827\n",
      "Energy expected   for eigenstate n= 7  is  241.805307827\n",
      "Energy calculated for eigenstate n= 8  is  315.827340835\n",
      "Energy expected   for eigenstate n= 8  is  315.827340835\n",
      "Energy calculated for eigenstate n= 9  is  399.718978244\n",
      "Energy expected   for eigenstate n= 9  is  399.718978244\n",
      "Energy calculated for eigenstate n= 10  is  493.480220054\n",
      "Energy expected   for eigenstate n= 10  is  493.480220054\n"
     ]
    }
   ],
   "source": [
    "# Define constants - here I use atomic units\n",
    "hbar = 1.0\n",
    "hbar_squared = hbar**2\n",
    "m_e = 1.0 # Mass of the electron\n",
    "two_m_e = 2.0*m_e\n",
    "# These are not needed but are easier\n",
    "pi_squared = np.pi**2\n",
    "square_well_width_squared = square_well_width**2\n",
    "\n",
    "energy = np.zeros(number_of_basis_functions)\n",
    "for m in range(number_of_basis_functions):\n",
    "    #range starts from 0 not 1 so need to add one to get n.\n",
    "    n = m+1\n",
    "    H_acting_on_phi = -(hbar_squared/two_m_e)*second_derivative_basis_functions[m]\n",
    "    energy[m] = integrate_functions(basis_functions[m],H_acting_on_phi,number_of_x_points,x_spacing)\n",
    "    print \"Energy calculated for eigenstate n=\",n,\" is \",energy[m]\n",
    "    print \"Energy expected   for eigenstate n=\",n,\" is \",hbar_squared*n*n*pi_squared/(two_m_e*square_well_width_squared)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These numbers agree extremely well! We would find worse agreement if we used fewer points in our x array (you might like to try this).\n",
    "\n",
    "### Matrix representation\n",
    "\n",
    "What about the rest of the Hamiltonian matrix, i.e. the numbers $\\langle \\phi_n \\vert \\hat{H} \\vert \\phi_m \\rangle$ ? We'll calculate these using two loops, and output the results.  I will use two approaches: a formatted print written in the loop, and the internal numpy printing for arrays (though this requires some changes to the options to make it look nice)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output of the Hamiltonian matrix elements as we calculate them\n",
      "   4.935   -0.000    0.000   -0.000    0.000    0.000    0.000   -0.000    0.000   -0.000\n",
      "  -0.000   19.739   -0.000   -0.000    0.000   -0.000    0.000   -0.000   -0.000   -0.000\n",
      "   0.000    0.000   44.413    0.000   -0.000    0.000   -0.000   -0.000    0.000    0.000\n",
      "  -0.000   -0.000    0.000   78.957    0.000   -0.000   -0.000   -0.000   -0.000   -0.000\n",
      "  -0.000   -0.000    0.000   -0.000  123.370    0.000    0.000    0.000   -0.000   -0.000\n",
      "   0.000    0.000    0.000    0.000    0.000  177.653   -0.000   -0.000    0.000    0.000\n",
      "   0.000    0.000   -0.000   -0.000    0.000    0.000  241.805   -0.000    0.000   -0.000\n",
      "  -0.000    0.000   -0.000    0.000    0.000   -0.000   -0.000  315.827   -0.000    0.000\n",
      "   0.000   -0.000    0.000   -0.000   -0.000    0.000    0.000   -0.000  399.719    0.000\n",
      "  -0.000   -0.000    0.000   -0.000   -0.000    0.000   -0.000    0.000    0.000  493.480\n",
      "\n",
      "Now the numpy output with formatting\n",
      "\n",
      "[[   4.935   -0.       0.      -0.      -0.       0.       0.      -0.       0.      -0.   ]\n",
      " [  -0.      19.739    0.      -0.      -0.       0.       0.       0.      -0.      -0.   ]\n",
      " [   0.      -0.      44.413    0.       0.       0.      -0.      -0.       0.       0.   ]\n",
      " [  -0.      -0.       0.      78.957   -0.       0.      -0.       0.      -0.      -0.   ]\n",
      " [   0.       0.      -0.       0.     123.37     0.       0.       0.      -0.      -0.   ]\n",
      " [   0.      -0.       0.      -0.       0.     177.653    0.      -0.       0.       0.   ]\n",
      " [   0.       0.      -0.      -0.       0.      -0.     241.805   -0.       0.      -0.   ]\n",
      " [  -0.      -0.      -0.      -0.       0.      -0.      -0.     315.827   -0.       0.   ]\n",
      " [   0.      -0.       0.      -0.      -0.       0.       0.      -0.     399.719    0.   ]\n",
      " [  -0.      -0.       0.      -0.      -0.       0.      -0.       0.       0.     493.48 ]]\n"
     ]
    }
   ],
   "source": [
    "# Create a square array to store the matrix elements\n",
    "H_matrix = np.eye(number_of_basis_functions)\n",
    "\n",
    "print \"Output of the Hamiltonian matrix elements as we calculate them\"\n",
    "# Loop over basis functions phi_n (the ket in the matrix element)\n",
    "for n in range(number_of_basis_functions):\n",
    "    # Act with H on phi_n and store in H_acting_on_phi_n\n",
    "    # This is another ket\n",
    "    H_acting_on_phi_n = -(hbar_squared/two_m_e)*second_derivative_basis_functions[n]\n",
    "    # Loop over basis functions phi_m (the bra in the matrix element)\n",
    "    for m in range(number_of_basis_functions):\n",
    "        # Create matrix element by integrating\n",
    "        H_matrix_element_mn = integrate_functions(basis_functions[m],H_acting_on_phi_n,number_of_x_points,x_spacing)\n",
    "        H_matrix[m,n] = H_matrix_element_mn\n",
    "        # The comma at the end prints without a new line; the %8.3f formats the number\n",
    "        print \"%8.3f\" % H_matrix_element_mn,\n",
    "    # This print puts in a new line when we have finished looping over m\n",
    "    print\n",
    "\n",
    "print\n",
    "print \"Now the numpy output with formatting\"\n",
    "print\n",
    "# Another way to output is this, though it's commented out to avoid clutter\n",
    "np.set_printoptions(precision=3,linewidth=100,suppress=True)    \n",
    "print H_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should not surprise you that the matrix is diagonal: because the basis vectors are eigenvectors, an operator acts on a basis vector and returns that basis vector multiplied by the relevant eigenvalue.  The orthogonality of the basis then ensures that only the diagonal elements are non-zero.\n",
    "\n",
    "Of course, we can choose other basis vectors than the eigenvectors.  In these cases, the matrix will *not* be diagonal.  We will look at this in a subsequent notebook.\n",
    "\n",
    "Once we have the Hamiltonian (or any observable matrix) we can find the eigenvalues and eigenvectors from the matrix.  There are many approaches to this, but we will rely on the built-in `numpy` routines in `numpy.linalg`: `eigh` for a Hermitian matrix (or just `eig` for a non-Hermitian matrix)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Overall procedure\n",
    "\n",
    "In general, when using the matrix approach to quantum mechanics, we will follow this pattern:\n",
    "\n",
    "* Choose a basis for the problem, ${\\vert\\phi_{m}\\rangle}$\n",
    "* Create the Hamiltonian matrix, $\\langle \\phi_{m} \\vert \\hat{H} \\vert \\phi_{n} \\rangle$\n",
    "* Find the eigenvalues and eigenvectors of the Hamiltonian\n",
    "* Or just evaluate the energy of a particular state, $\\vert \\psi \\rangle$\n",
    "\n",
    "There are many numerical libraries which allow you to solve for the eigenvalues and eigenvectors of a matrix.  We will introduce the `numpy` version of these, and show some examples, in the next notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A general state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is helpful to remember that we can also explore what happens to a general state, written as:\n",
    "\n",
    "$$\\vert \\psi \\rangle = \\sum_n c_n \\vert \\phi_n \\rangle$$\n",
    "\n",
    "We define an array of coefficients, and then build $\\vert\\psi\\rangle$ and its second derivative, and evaluate the energy in two ways: first, directly, as the expectation value of the Hamiltonian operator in this state ($\\langle \\psi \\vert \\hat{H} \\vert \\psi \\rangle$); second, by summing over the eigenvalues, weighted by the square modulus of the coefficients ($E = \\sum_n \\vert c_n \\vert^2 E_n$).  Of course, these two should be exactly the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy numerically is:  43.5907527715\n",
      "Energy analytically is:  43.5907527715\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu8lXP2wPHPOl3dcpskSUZCUWQqRTMdl1RIF4nIJWka\nRjWjMcKM4jejaSYzIUkcjEsiKSGU6mjSZaSr7pnShcJQKpxS398f68SZnPt+9v4+l/V+vfarvc95\nzvOsntc5e+3vZX2/4pzDGGNMMmX5DsAYY4w/lgSMMSbBLAkYY0yCWRIwxpgEsyRgjDEJZknAGGMS\nLOUkICLHisg0EVkqIktEpG8Rxz0oIqtFZKGInJHqdY0xxqSuYgDn+A641Tm3UEQOBt4XkcnOuRX7\nDhCRdkBd51w9ETkLGAk0D+DaxhhjUpByS8A5t9k5tzD/+Q5gOVBrv8M6AE/nHzMXOFREaqR6bWOM\nMakJdExARI4HzgDm7vetWsCGAq838eNEYYwxJsMCSwL5XUEvAf3yWwTGGGNCLogxAUSkIpoAnnHO\nvVLIIZuA2gVeH5v/tcLOZYsZGWNMGTnnpDw/F1RL4AlgmXPugSK+PxG4FkBEmgNbnXNbijqZc84e\nzjFw4EDvMYThYffB7oXdi+IfqUi5JSAi5wBXA0tEZAHggDuBOvp+7kY55yaJyEUisgbYCfRI9bph\n9v77sHYt5OXp4+ij4cILoWIg7S5jjAlOym9Lzrl3gQqlOO6WVK8VZnv2wCuvwNChsHkzNG4MVatC\nlSqwYgXceCNcfz307Al16/qO1hhjlH02DcDChdC1KxxxBPzud9CpE1TYLy0uWwY5OdC8Odx5J/zm\nNyAl9OBlZ2enLeYosfvwA7sXP7B7EQxJtT8paCLiwhZTcd56C665Bh56SBNBSW/s69ZB585Qvz48\n9hgceGBGwjTGxJiI4DwPDCdSTg5cdx28/DJccUXJCQDg+OPh3Xe1pXD22bBhQ4k/YowxaWMtgXIa\nMQLuvx8mTYKTTy77zzsHgwfDmDEwcyZUqxZ8jMaYZEilJWBJoBxmz4YOHWDOHDjhhPKfxzn41a9g\n40YdVLbZQ8aY8rDuoAz67DPt+snJSS0BgHYfDR8Ou3bpgLIxxmSaJYEy2LMHrroKuneH9u2DOWel\nSjB2rA4wP/poMOc0xpjSsu6gMhg0CGbMgMmTg++6Wb0aWrSAuXOtjsAYUzY2JpABK1ZAy5awZAnU\nrJmea+wbaH777dLNNDLGGLAxgbRzTou77rorfQkAoF8/+OoreOKJ9F3DGGMKspZAKbz6Ktx+Oyxa\npH346bR4MZx/vl7rmGPSey1jTDxYd1AaffstnHYaPPIItG6dmWv+4Q+6zMTLL2fmesaYaLPuoDT6\nxz80CWQqAYAmgcWLYerUzF3TGJNM1hIoxpYtcOqp8O9/p14TUFbPPw/DhmlBmg0SG2OKYy2BNBk2\nDK68MvMJALQgLS8PJkzI/LWNMclhLYEibN0KJ54I8+bpom8+TJoEt92mXUP7L01tjDH7WEsgDUaM\ngIsv9pcAANq10z0Knn3WXwzGmHizlkAhvv5au4CmTYMGDbyGwsyZukzFypW6S5kxxuzPWgIBy8nR\nJRx8JwDQKuVTT9WYjDEmaNYS2M/u3ToWMHYsNGvmLYz/8e67unnNypU2NmCM+TFrCQRozBioVy88\nCQB0B7Lq1XXPAWOMCZIlgf08/LCu4RMmIrrfwNChviMxxsSNJYECFi2CTZt0Vk7YdOyoxWuzZvmO\nxBgTJ5YEChg1Cm68MZzbPFaoALfeaq0BY0ywbGA4386dULu2FmYde2zGL18qO3dq3cKsWTpuYYwx\nYAPDgRgzRqdjhjUBABx0EPTurYvaGWNMEKwlkK9ZMxg4UKuEw+yTT7R+Yf16OOQQ39EYY8LAWgIp\nWrAANm+Gtm19R1KymjXh3HNh9GjfkRhj4iCQJCAiOSKyRUQWF/H9ViKyVUTm5z/+EMR1g/Loozog\nHJVCrN69NeaQNeKMMREUSHeQiLQEdgBPO+caFfL9VkB/59ylpThXRruDvv1Wt3FctEgHhqNg716t\nan7hBWja1Hc0xhjfvHcHOedmAl+WcFgot0Z54w04/fToJACArCzo1UtbA8YYk4pMjgm0EJGFIvK6\niIRgaTY1ejRcdZXvKMquRw8YNw62bfMdiTEmyjKVBN4HjnPOnQEMB0KxX9ZXX8HkyXDZZb4jKbuj\nj4YLLoDnnvMdiTEmyjJSG+uc21Hg+RsiMkJEjnDOfVHY8YMGDfr+eXZ2NtnZ2WmJa/x4yM7WjVui\nqHdv6N8fbrrJ9iE2Jklyc3PJzc0N5FyB1QmIyPHAq865hoV8r4Zzbkv+82bAi86544s4T8YGhtu0\ngRtu0P18o2jfAPGLL0KTJr6jMcb4ksrAcCAtAREZDWQDR4rIemAgUBlwzrlRQBcRuQnYDXwDeH/b\n3bIF5s7V1kBUZWXBtdfCM89YEjDGlE9iK4aHD9ck8Mwzab9UWq1ZA+ecAxs3QqVKvqMxxvjgfYpo\nFEV1VtD+TjxRH2+95TsSY0wUJTIJrFsHq1fr7Jo42NclZIwxZZXIJDB+PHToEJ/uk65d4c03YetW\n35EYY6ImkUlgwgTdqSsuDj8cWreGl17yHYkxJmoSlwQ++0zXCYpLV9A+11wDTz/tOwpjTNQkLglM\nnAgXXghVq/qOJFjt2sHy5bB2re9IjDFRkrgkELeuoH0qV9aiN1tGwhhTFomqE9i+HWrVgg0b4NBD\n03IJr2bO1CUklizxHYkxJpOsTqCU3noLWrSIZwIAOPtsnSG0bJnvSIwxUZGoJDB+PHTq5DuK9MnK\ngssv17WEjDGmNBLTHbRrly6//MEHupNYXM2Zo3sNLFtmK4sakxTWHVQKublw8snxTgAAZ50F33xj\n4wLGmNJJTBKYOFGrhONORCuIX3jBdyTGmChIRBJwTvcSvvhi35FkxhVXaBIIWU+fMSaEEpEEVq2C\nvDw47TTfkWTGmWfqv/Pn+43DGBN+Gdle0rdJk+Cii5IzUCryQ2vgZz/zHY1Jor17dRLGihXw+ee6\nXMvXX0ONGjouV6uWflg56CDfkZpEzA668EK4+eZ4VgoXZdEi/f/+5z/JSX7Gr2+/hbFjtet16lSt\nx2nUCKpX18cBB+iOfh9/DOvXw9KluiNe69bQuTOccorv/0F0pTI7KPZJYMcOqFlTf/EOOSSw04ae\nc1Cvnv5RNm7sOxoTZ1u3wsiR8OCD+qZ/+eW6QGOdOsX/3I4dMGMGTJ6srdaGDaFvX221ZyWiozo4\nNkW0GNOmQbNmyUoAoJ/+O3eGl1/2HYmJqz17YNgwqFtX61LefFMfPXuWnAAADj5Y3/CHDdONnq69\nFu65B+rXh1desYkNmRL7lsBNN+n2i/37B3bKyJgzR/8gly71HYmJmxUr4IYbdGOmxx/XVmcQnIMp\nU+C3v9Wxg3/8IzkTOlJhLYEiOKeDwu3a+Y7Ej2bNtKm+YoXvSExcOAfDh0PLlrpH9/TpwSUA0Bbs\nhRfCwoVa13PeefDHP8Lu3cFdw/yvWCeBfUsn1K/vOxI/srJ0rSTrEjJB+O476NMHHnkE5s6FW25J\nX999pUp6/sWLYd48TTqrV6fnWkkX6yTwxhvaCkjy7BgbFzBB2L5dP5mvWgWzZuk4QCYcfbS25rt3\n11Vyn3kmM9dNklgngX31AUn2i1/ARx/pw5jy+OILaNVK5/a//nrml2IX0RbItGk6cHzbbToobYIR\n2ySwcye89x6ce67vSPyqWBEuvVSX0TamrLZvh7ZttW/+0Ue1m8aXhg21G+r99/V3+quv/MUSJ7FN\nAv/6l1YkHnyw70j869wZxo3zHYWJmm++gfbttc7kb38LR7fqkUfq5lB16mj30KZNviOKvtgmgSlT\ntGDF6H1YsgQ2b/YdiYmKXbugSxftAhoxIhwJYJ9KlTSm667TAeM1a3xHFG2xTQJvv63l6AaqVIE2\nbeC113xHYqKiXz9943/qKahQwXc0hbvtNrjzTh2vWLzYdzTRFcsksGWLDoQ2aeI7kvDo0EGrMI0p\nyeOP6/z/0aP9jgGURq9eWlDWujX8+9++o4mmQJKAiOSIyBYRKTIfi8iDIrJaRBaKyBlBXLcoU6dC\ndrYOihrVrh28844OmBtTlDlz4I47YMIEqFbNdzSl07Ur5OTo+IUtn152QbUEngTaFPVNEWkH1HXO\n1QN6AyMDum6hrCvoxw4/HJo21cW6jCnM5s26+FtOTvRW9LzkEl3E7qKLrGuorAJJAs65mcCXxRzS\nAXg6/9i5wKEiUiOIa/84FhsULop1CZmi7N2ry0DccINOv4yiTp10JdM2bWy9rLLI1JhALWBDgdeb\n8r8WuFWrdEDrpJPScfZo69BBi32s0Mbsb9gwnRF0992+I0lN1646nbVtWyuQLK1Q9poPGjTo++fZ\n2dlkZ2eX+mf3tQLCNKUtLOrU0Sl/s2bBz3/uOxoTFkuXwuDBWogV1plAZdG9u1Y5t2kDM2fCT37i\nO6Lg5ebmkpubG8i5AltKWkTqAK865xoV8r2RwHTn3Av5r1cArZxzWwo5NqWlpDt21K0Vu3Ur9yli\nbeBAHRweOtR3JCYMdu2Cs87Sxdp69vQdTbDuukvHB6dOjX/RaFiWkpb8R2EmAtcCiEhzYGthCSBV\n330Hublw/vlBnzk+9o0LhGwbCePJoEFQu7aOBcTNn/6kS0106WJLURcnqCmio4FZwEkisl5EeohI\nbxH5JYBzbhKwVkTWAI8CNwdx3f3Nm6ddHkcdlY6zx0PjxpCXB8uX+47E+LZggc4EeuyxeHafiuiM\noYoV4de/tg8+RYnVzmKDB2uh2LBhAQcVM7fcomMDd9zhOxLjy5490KKF7rzXo4fvaNJrxw4dA+vW\nDX7/e9/RpEdYuoO8mz7dVg0tjfbtbQmJpBs1CqpW1fV34u7gg+HVV+Ghh+Cll3xHEz6xaQns2qUr\nDK5fr4VRpmh5eVCjhu7UVL2672hMpm3erH3lublw6qm+o8mcBQt068rXXtPB8DixlgC6d0C9epYA\nSqNKFR08nzTJdyTGh/794cYbk5UAQMfDcnLgsstsCeqCYpMErCuobNq31yaySZbp0+Hdd3Xz9iS6\n9FIdJO7YUfdLMDFKArm5umicKZ2LLtI51Hl5viMxmbJnD/z2t1ojcuCBvqPxZ8AAXVGgZ0+bMQQx\nSQJ5eVrt+Itf+I4kOo46CurXhxkzfEdiMuXpp3WQ9LLLfEfil4gul716NfzlL76j8S8WSWDuXF31\nMNMbYEeddQklx86d2gV0//3xrAkoqwMO0OWyhw+HN9/0HY1fsUgC1hVUPvuSgDWJ42/oUJ0rH7dZ\nMamoVQteeEGnySZ5i8pYJAEbFC6f007TBGDL7sbbxx/rEsuDB/uOJHxattT1tDp10qKyJIp8ncC3\n3+oqgR9/HJ2dkMKkTx845hirHo6zG2/UGpohQ3xHEk7O6T3avl1bBlHsLkt0ncCcOTrf2RJA+di4\nQLytXq193wMG+I4kvETg4Ydh3TodM0mayCcBGw9ITatW2h30+ee+IzHpcM898JvfWBFlSapW1SUl\nhg7V7uUkiXwSeOcdSwKp2Fc9/MYbviMxQVu6VDdZ6tfPdyTRcNxx8Oyzus3mxo2+o8mcSCeBXbt0\nuYizz/YdSbRdcoktKBdHAwfCbbfBIYf4jiQ6LrgA+vaFyy/X95ckiPTA8KxZuizy/PlpDirmNm+G\nBg10Ge5KlXxHY4Iwf74m9zVrkl0dXB5790LnzrrZzkMP+Y6mdBI7MDxjhlUJB+Hoo6FuXV1TxsTD\n3XfDnXdaAiiPrCx46intIh0zxnc06RfpJPCvf1kSCIp1CcXHvHmwaBH06uU7kug67DAdKO7TB5Yt\n8x1NekU2CezZo59cW7b0HUk8WBKIj/vu07GAKlV8RxJtZ5wBf/2rrrW0fbvvaNInsklg8WKoWdP2\nEw5K48bw1VfJLp+Pgw8+0LGyG2/0HUk89OgB55yjraqQDZ8GJrJJwMYDgpWVBRdfDK+/7jsSk4rB\ng7UuwMYCgvPQQ7BypRaUxVFkk4CNBwTv4outSyjK1qyBt96Cm2/2HUm8HHCAjg/ce6+uWBw3kZwi\n6pzukTtvnhZ4mGDs2KFdbJs22TIcUdSrl870+r//8x1JPE2YoIV38+frWkxhkrgpoitXanPXEkCw\nDj5Y+z+nTPEdiSmrjRth3DirDk6njh2ha1fo3l1rCeIikknAxgPSx2YJRdM//gHXX68r6pr0ue8+\nbTHfd5/vSIITye6ga67Rhc9sBkTw1q6F5s3hk090sNiE39atcMIJWhtQu7bvaOJv0yZo0kTXGTr/\nfN/RqMR1B82YobskmeD99Kf6aXLePN+RmNIaNQouusgSQKbUqqUJoHt3TQhRF7kksGEDfPMNnHSS\n70jiy7qEomPXLnjgAS0OM5lz/vnw61/DlVfC7t2+o0lN5JLAu+/q4GUUd/+JCksC0fH887qp0umn\n+44kee68UydTRH1XvkCSgIi0FZEVIrJKRG4v5PutRGSriMzPf/yhvNeaOdOWiki3Fi3go4/i0dSN\nM+d0E5Tf/c53JMmUlaXdQmPHwvjxvqMpv5STgIhkAcOBNsCpQDcROaWQQ2c4587Mf/ypvNezJJB+\nFStCmzZWPRx2b70FFSpA69a+I0muI4/UJNC7d3SXXAmiJdAMWO2c+8g5txsYA3Qo5LiUO3C2bdMb\n3bhxqmcyJbnkEksCYTd0KPTvb12jvjVrphv4dOmi45VRE0QSqAVsKPB6Y/7X9tdCRBaKyOsi0qA8\nF5ozR6dmVa5cnp82ZdG2re61GsVf6iRYsgSWL4crrvAdiQFdqqN+fd3kKmoqZug67wPHOee+FpF2\nwASgyPk9gwYN+v55dnY22fmbCFtXUOYccYQONubmQrt2vqMx+3vwQbjpJvtAFBYi8Nhj2irIyYGe\nPdN7vdzcXHJzcwM5V8rFYiLSHBjknGub/3oA4JxzQ4r5mbXAz5xzXxTyvSKLxc49F26/XT+lmvQb\nMkQHiEeM8B2JKeizz3SK9KpVUL2672hMQcuX62oGkydnttvad7HYe8CJIlJHRCoDVwIT9wuwRoHn\nzdDk86MEUJzdu7WAqUWLACI2pbJvqmjIisoTb9Qo3ejEEkD41K8Pw4fr+MCXX/qOpnRSTgLOuT3A\nLcBkYCkwxjm3XER6i8gv8w/rIiIfiMgCYBhQ5p7MBQu0NP7QQ1ON2JRWgwY6U2jJEt+RmH127dKW\nmS0UF15XXKEfoK67LhoLzUVm7aC//x0+/DC+GzuEVb9+unvbXXf5jsQAPPccPPEETJ3qOxJTnF27\nIDsb2rfPTDGZ7+6gjLBBYT+sejg8nINhw3TnMBNulSvDiy/qrmRvv+07muJFoiXgnG6WMW+eLZKV\nabt2aUtg1Srbz9m32bN1Bd1Vq2yF16jIzYVu3XRHsnTufxL7lsCaNVCliiUAHypXhgsugEmTfEdi\nhg/XRcssAURHdrYW9HXpAnl5vqMpXCR+nfYtGmf8aN8eXn3VdxTJtnmzJuIePXxHYsqqf39tBfTt\n6zuSwlkSMCVq1077NcP6SSYJRo3SWSeHHeY7ElNWIvDkk7oPSk6O72h+LBJJYNYsOPts31Ek11FH\n6XLF77zjO5Jk2rULRo7UriATTYccoiuNDhgQvg2bQp8EvvhCN5Jp1Mh3JMlmXUL+jB8PJ58MDRv6\njsSk4pRT4NFHtdDvs898R/OD0CeBOXOgaVMtWjL+XHKJJoGQTSZLhIceiubCZObHOnfW2UJXXgnf\nfec7GhX6JGDjAeFw2mnat2nVw5m1cKGu39ShsMXZTST9+c+6D0RYdiQLfRKw8YBwENE3ookTSz7W\nBOfhh+FXv7KWcJxUqKDbgo4bB2PG+I4m5MViu3frksYbNtisiDCYNk1XcX3vPd+RJMOXX+p6WStW\nQI0aJR9vomXRIq3BmTo19THP2BaLLVoExx9vCSAsfv5zXb/J9h7OjH/+U6fnWgKIp9NP130hOnXS\nCTC+hDoJ2HhAuFSqpG9KNkso/fbu1dVCbVpovHXrpknA50BxqJOAjQeEj40LZMbUqXDAAfb7nwR/\n+YvOuvM1UBzaJOCctQTCqG1bXdF1+3bfkcTbiBG6b61tIh9/FSvqAPG4cTB6dOavH9oksGGDDgyf\ncILvSExB1arp7m6TJ/uOJL42bNDq7Kuv9h2JyZQjj4QJE3T/jvnzM3vt0CaBfa0A+yQUPh06wCuv\n+I4ivh59FLp3h4MP9h2JyaRGjeCRR7Sg7NNPM3fd0CaBWbNsP+Gwat9eV7QMS8VjnOzaBY8/Djfd\n5DsS40OXLrpnRJcu+ruQCaFOAjYeEE61a0OdOtpaM8EaP173dq5f33ckxpd77tFp8ZnaRzqUSWDH\nDi2QOfNM35GYonTqpG9YJlgjRlgrIOmysuDZZ3VcaOTIDFwv/Zcou/fe00KKqlV9R2KKsi8JhKzg\nPNI++ABWr4aOHX1HYnyrVk3H3QYO1H0I0imUSWD2bJsfHXYNGuiWn5meyRBnI0dCr15alGdMvXra\nIrjiCli3Ln3XCWUSsEHh8BPR1sDLL/uOJB527NA54r16+Y7EhEnr1roRTYcO+juSDqFMArNnWxKI\ngs6dbVwgKM89p5uSH3us70hM2PTtC02awLXX6nIiQQtlEqhWDY45xncUpiRNm8K2bTqIb8rPuR8q\nhI3Zn4j+fmzZojOHghbKJGDjAdGQlaWDmNYaSM2sWfDNN3Deeb4jMWFVpYp2vT71FIwdG+y5LQmY\nlFiXUOr2tQKyQvnXaMKiRg1dWuLmm4OdkBHKTWXmz3c0buw7ElMau3fD0UfrNoi1a/uOJno+/VQ3\nkf/Pf+Dww31HY6LgpZfg1lvh3//Wvz0IwaYyItJWRFaIyCoRub2IYx4UkdUislBEzijufA0bBhGV\nyYRKlXQTemsNlE9ODlx2mSUAU3pdusANN2grPC8v9fOlnAREJAsYDrQBTgW6icgp+x3TDqjrnKsH\n9AaKrYOz/VSj5bLLdBlcUzZ79mhtgA0Im7K6+26dPPPLX6ZesBlES6AZsNo595FzbjcwBuiw3zEd\ngKcBnHNzgUNFxDbNi4kLL4TFi+GTT3xHEi2vv65/yLY8iimrrCzdfnTxYrj//hTPFUA8tYANBV5v\nzP9accdsKuQYE1FVq2qXkLUGysamhZpUHHSQLi3x97+ndp5QdrwMGjTo++fZ2dlkZ2d7i8WUTteu\n8Le/wS23+I4kGtas0RkeEyb4jsREUW5uLrm5uYDu+/3EE+U/V8qzg0SkOTDIOdc2//UAwDnnhhQ4\nZiQw3Tn3Qv7rFUAr59yWQs7nwjZjyZQsL09nKixdaoV+pdG/v459DRlS8rHGlMT37KD3gBNFpI6I\nVAauBPbfinwicC18nzS2FpYATHRVqaKbzViXUMl27tT+XFsy2oRByknAObcHuAWYDCwFxjjnlotI\nbxH5Zf4xk4C1IrIGeBSwntAYuvzy4KsZ42j0aN0w6fjjfUdiTEiLxcIWkymdvDyoWVPXxbcuocI5\nB2ecAUOH6gqRxgTBd3eQMYB2CdksoeLNnKnJ8vzzfUdijLIkYALVtSu8+KLvKMJr+HD49a9tnSAT\nHtYdZAKVl6ddQQsWwHHH+Y4mXDZt0iVR1q3T5dKNCYp1B5nQqFJFl5F4/nnfkYTPqFHQrZslABMu\n1hIwgXvnHejTR0vajcrLgzp1YNo03Z/ZmCBZS8CEys9/Dlu3wpIlviMJjzFjoFEjSwAmfCwJmMBl\nZcFVV+m+uUanhT7wAPTr5zsSY37MkoBJi6uv1qKodGyMHTUzZ8KOHbrGizFhY0nApEXDhnDYYfCv\nf/mOxL8HHoC+fW1aqAkn+7U0aXP11dYltG4dTJ8O11/vOxJjCmezg0zarF8PjRvDxx/r1NEkuu02\n7RJLdeMPY4pjs4NMKB13nHYLvfaa70j82LEDnnzS9lgw4WZJwKTVDTektuFFlD35JLRqBT/9qe9I\njCmadQeZtPr6azj2WK0ZqJWgDUW/+w7q1dPK6ebNfUdj4s66g0xoHXig7jPw1FO+I8msceM0+VkC\nMGFnScCkXc+e2iWUlJoB53S/5d/9znckxpTMkoBJu6ZNtUUwY4bvSDLjnXdg+3bdbtOYsLMkYNJO\nRFsDOTm+I8mMoUN1I3krDjNRYAPDJiM+/xxOPFGLpw47zHc06bNsGZx3nv4/q1b1HY1JChsYNqH3\nk5/AhRfGf5+Bv/5Vdw6zBGCiwloCJmPefhtuvRUWLdIuorhZuxaaNIEPP4x3a8eEj7UETCScfz7s\n3h3fAeIhQ+BXv7IEYKLFWgImox5+WBdUe+kl35EEa+NG3TRm5UqoXt13NCZpUmkJWBIwGbV9Oxx/\nPCxcCLVr+44mOL/5DVSoYAvFGT8sCZhI6dcPDjoI7rvPdyTB2LIF6teHDz6AY47xHY1JIksCJlJW\nrYKWLeGjj+CAA3xHk7oBA7SF8/DDviMxSWVJwEROu3bQtSv06OE7ktR8+qm2AubPhzp1fEdjkspm\nB5nI6dMHHnpI19mJsvvu0x3ULAGYqEqpJSAihwMvAHWAdUBX59y2Qo5bB2wD9gK7nXPNijmntQQS\nYO9eOO003X+3dWvf0ZTPvp3Tli2DGjV8R2OSzGdLYADwtnPuZGAacEcRx+0Fsp1zjYtLACY5srLg\nzjvhz3/2HUn53Xuv1gVYAjBRlmpLYAXQyjm3RUSOBnKdc6cUctxaoIlz7r+lOKe1BBLiu+/g5JPh\nn//UgeIoWblSY161Cg4/3Hc0Jul8tgSOcs5tAXDObQaOKuI4B0wRkfdEpFeK1zQxUbEi3H57NFsD\nd9+tK4VaAjBRV2JLQESmAAUbvIK+qf8BeMo5d0SBY//rnDuykHPUdM59IiLVgSnALc65mUVcz1oC\nCZKXB3XrwiuvwM9+5jua0pk3Dy69FFav1noHY3xLpSVQsaQDnHNFDtuJyBYRqVGgO+jTIs7xSf6/\nn4nIeKAZUGgSABg0aND3z7Ozs8nOzi4pTBNRVaroDlz33adbMoadc9C3L/zpT5YAjD+5ubnk5uYG\ncq5UxwT5/iYAAAAHgElEQVSGAF8454aIyO3A4c65AfsdcyCQ5ZzbISIHAZOBe5xzk4s4p7UEEmbn\nTjjhBJg2DU491Xc0xXv2WZ3RNHeubRpjwsNbsZiIHAG8CNQGPkKniG4VkZrAY865S0Tkp8B4tAup\nIvCcc+4vxZzTkkAC/f3vurDcq6/6jqRo27drYdjYsdCihe9ojPmBVQybyMvLgwYN4LHHdGeuMLrj\nDti0CZ5+2nckxvwvSwImFsaO1bGBefN0Rc4wWbMGmjeHxYttkTgTPrZshImFLl10QblnnvEdyf9y\nTreMvO02SwAmfqwlYEJl9my4/HItxgrL7JvHH4eRI2HOHK1tMCZsrDvIxMoVV+gAbIGZwt6sX6/1\nC9On61pHxoSRJQETK/veeKdNg4YN/cXhHLRpA61awV13+YvDmJLYmICJleOOg8GD4frrdWN6X3Jy\n4IsvdGkLY+LKWgImlJzTjWfOOQf++MfMX3/ZMm0BWDeQiQLrDjKxtHGjrtf/9ttw+umZu+62bdC0\nqS51ff31mbuuMeVlScDE1pNP6jINs2dnZj/ivXuhY0eoXdv2DDbRYUnAxJZz0L27jg2MGZP+9Xru\nvRcmT9ZB6cqV03stY4JiA8MmtkR0gHbTJl3DP52efx5GjdLKZUsAJims9MWEXtWqMGECnHUW1KsH\n110X/DXGjIFbb4UpU6BmzeDPb0xYWRIwkVC9Orz2GmRn6/OLLgru3GPGwG9/qwnAZgKZpLHuIBMZ\nDRroDmQ9e8KIEcGc84knLAGYZLOBYRM5H34Il1wCbdvC0KHlW3F02za4+WZYsABeekkTjDFRZQPD\nJlHq1oVZs3RZ59at9Y28LGbN0vqDatV02WpLACbJLAmYSDr8cHjzTV1++uKLddG5VauKPn7PHh1c\nvuACuOwy3cnskUfgwAMzF7MxYWTdQSbydu7UgrL774cjjoAmTbTit2JF2LBBK49nzdK9APr00SRQ\npYrvqI0JjhWLGYN+2l+xQrt45s3T6t/atfXRsCE0auQ7QmPSw5KAMcYkmA0MG2OMKRdLAsYYk2CW\nBIwxJsEsCRhjTIJZEjDGmASzJGCMMQlmScAYYxLMkoAxxiRYSklARLqIyAciskdEzizmuLYiskJE\nVonI7alc0xhjTHBSbQksAToB7xR1gIhkAcOBNsCpQDcROSXF6yZCbm6u7xBCwe7DD+xe/MDuRTBS\nSgLOuZXOudVAceXKzYDVzrmPnHO7gTFAh1SumxT2S67sPvzA7sUP7F4EIxNjArWADQVeb8z/mjHG\nGM9K3GNYRKYANQp+CXDAXc65V9MVmDHGmPQLZBVREZkO9HfOzS/ke82BQc65tvmvBwDOOTekiHPZ\nEqLGGFNG5V1FtMSWQBkUFcB7wIkiUgf4BLgS6FbUScr7HzHGGFN2qU4R7SgiG4DmwGsi8kb+12uK\nyGsAzrk9wC3AZGApMMY5tzy1sI0xxgQhdJvKGGOMyRwvFcOlKR4TkQdFZLWILBSRMzIdY6aUdC9E\n5CoRWZT/mCkiDX3EmQmlLSoUkaYisltEOmcyvkwq5d9ItogsyC/YnJ7pGDOlFH8j1URkYv57xRIR\nud5DmBkhIjkiskVEFhdzTNneO51zGX2giWcNUAeoBCwETtnvmHbA6/nPzwLmZDrOEN2L5sCh+c/b\nJvleFDhuKvAa0Nl33B5/Lw5Fu1dr5b/+ie+4Pd6LO4DB++4D8F+gou/Y03Q/WgJnAIuL+H6Z3zt9\ntARKUzzWAXgawDk3FzhURGoQPyXeC+fcHOfctvyXc4hvjUVpiwr7AC8Bn2YyuAwrzb24ChjnnNsE\n4Jz7PMMxZkpp7oUDDsl/fgjwX+fcdxmMMWOcczOBL4s5pMzvnT6SQGmKx/Y/ZlMhx8RBWQvpbgTe\nSGtE/pR4L0TkGKCjc+4Riq9Sj7rS/F6cBBwhItNF5D0RuSZj0WVWae7FcKCBiHwMLAL6ZSi2MCrz\ne2eQU0RNGonIuUAPtDmYVMOAgn3CcU4EJakInAmcBxwEzBaR2c65NX7D8qINsMA5d56I1AWmiEgj\n59wO34FFgY8ksAk4rsDrY/O/tv8xtUs4Jg5Kcy8QkUbAKKCtc664pmCUleZeNAHGiIigfb/tRGS3\nc25ihmLMlNLci43A5865b4FvRWQGcDrafx4npbkXPYDBAM65D0VkLXAKMC8jEYZLmd87fXQHfV88\nJiKV0eKx/f+IJwLXwvcVx1udc1syG2ZGlHgvROQ4YBxwjXPuQw8xZkqJ98I5d0L+46fouMDNMUwA\nULq/kVeAliJSQUQORAcB41h/U5p78RFwAUB+//dJwH8yGmVmCUW3gsv83pnxloBzbo+I7CseywJy\nnHPLRaS3ftuNcs5NEpGLRGQNsBPN9LFTmnsB/BE4AhiR/wl4t3Oumb+o06OU9+J/fiTjQWZIKf9G\nVojIW8BiYA8wyjm3zGPYaVHK34s/AU8VmDb5e+fcF55CTisRGQ1kA0eKyHpgIFCZFN47rVjMGGMS\nzLaXNMaYBLMkYIwxCWZJwBhjEsySgDHGJJglAWOMSTBLAsYYk2CWBIwxJsEsCRhjTIL9P0LQwuCN\nzQDqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x102f1d250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define an array to hold the coefficients of our state\n",
    "coefficient_array = np.zeros(number_of_basis_functions)\n",
    "coefficient_array[0] = 1.0\n",
    "coefficient_array[2] = 2.0\n",
    "coefficient_array[3] = 1.0\n",
    "# Prepare to make psi and its second derivative\n",
    "psi = np.zeros(number_of_x_points)\n",
    "second_derivative_psi = np.zeros(number_of_x_points)\n",
    "for i in range(number_of_basis_functions):\n",
    "    psi = psi+coefficient_array[i]*basis_functions[i]\n",
    "    second_derivative_psi = second_derivative_psi+coefficient_array[i]*second_derivative_basis_functions[i]\n",
    "# Now we will normalise psi and apply the same factor to d2psi\n",
    "integral = integrate_functions(psi,psi,number_of_x_points,x_spacing)\n",
    "psi = psi/np.sqrt(integral)\n",
    "second_derivative_psi = second_derivative_psi/np.sqrt(integral)\n",
    "# Also apply the normalisation to the coefficient array\n",
    "# We could also have normalised psi by ensuring that the square modulus of the\n",
    "# coefficient array was one\n",
    "coefficient_array = coefficient_array/np.sqrt(integral)\n",
    "# Plot psi\n",
    "pl.plot(x,psi)\n",
    "# Evaluate the energy from <psi|H|psi> - with no potential we just have\n",
    "# the kinetic energy term.  Find the expectation value explicitly\n",
    "H_acting_on_psi = -(hbar_squared/two_m_e)*second_derivative_psi\n",
    "print \"Energy numerically is: \",integrate_functions(psi,H_acting_on_psi,number_of_x_points,x_spacing)\n",
    "# Find the expectation value of the Hamiltonian (the energy) \n",
    "# by summing over energies weighted by the square modulus of the coefficients\n",
    "energy_psi = 0.0\n",
    "for i in range(number_of_basis_functions):\n",
    "    energy_psi = energy_psi + coefficient_array[i]*coefficient_array[i]*energy[i]\n",
    "print \"Energy analytically is: \",energy_psi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we see that the energies agree perfectly.  Notice that the numerical approach, of evaluating the expectation value of the Hamiltonian through integration on a grid, is completely general, and can be extended to situations where the potential changes (though in that case we will need to look at how many basis functions we should include)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}