{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Variational Method\n",
    "\n",
    "The variational method is a simple approach to finding the energy of the ground state of a system.  We write a *trial* wavevector which depends on one or more parameters, and then vary the parameters to find the lowest energy.  We will show in lectures that this always gives an *upper bound* to the energy of the system; we define the energy for the system in the trial state as:\n",
    "\n",
    "$$\n",
    "E(\\lambda) = \\frac{\\langle \\psi(\\lambda) \\vert \\hat{H} \\vert \\psi(\\lambda) \\rangle}{\\langle \\psi(\\lambda)\\vert \\psi(\\lambda)\\rangle}\n",
    "$$\n",
    "\n",
    "where $\\lambda$ is the parameter.\n",
    "\n",
    "We will start with a simple example: finding the optimum value of the exponent, $\\alpha$, for the ground state of the quantum harmonic oscillator.  As we know the exact answer, we can check that the approach is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "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 x-axis from -width to +width\n",
    "# This makes the QHO finite in width, which is an approximation\n",
    "# We will need to be careful not to make omega too small\n",
    "width = 10.0\n",
    "num_x_points = 1001\n",
    "x = np.linspace(-width,width,num_x_points)\n",
    "dx = 2.0*width/(num_x_points - 1)\n",
    "\n",
    "# Integrate two functions over the width of the well\n",
    "def integrate_functions(f1,f2,size_x,dx):\n",
    "    \"\"\"Integrate two functions over defined x range\"\"\"\n",
    "    sum = 0.0\n",
    "    for i in range(size_x):\n",
    "        sum = sum + f1[i]*f2[i]\n",
    "    sum = sum*dx\n",
    "    return sum\n",
    "\n",
    "# Define a function to act on a basis function with the potential\n",
    "def add_pot_on_basis(Hphi,V,phi):\n",
    "    for i in range(V.size):\n",
    "        Hphi[i] = Hphi[i] + V[i]*phi[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Define the potential in the square well\n",
    "def square_well_potential(x,V,a):\n",
    "    \"\"\"Potential for a particle in a square well, expecting two arrays: x, V(x), and potential height, a\"\"\"\n",
    "    for i in range(x.size):\n",
    "        V[i] = 0.5*a*x[i]*x[i]\n",
    "    # If necessary, plot to ensure that we know what we're getting\n",
    "    #pl.plot(x,V)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's define the trial function to use, which will be a gaussian:\n",
    "$$f(x) = A exp(-0.5\\alpha^2 x^2)$$\n",
    "\n",
    "We'll also need its second derivative for the energy, which we can find analytically in this case:\n",
    "$$df/dx = -x\\alpha^2 A exp(-0.5\\alpha x^2)$$\n",
    "$$d2f/dx2 = -\\alpha^2 A exp(-0.5\\alpha x^2) + \\alpha^4 x^2 A exp(-0.5\\alpha x^2) = \\alpha^2 A (\\alpha^2 x^2 -1) exp(-0.5\\alpha x^2)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "From theory, E min, alpha:  1.0 1.41421356237\n",
      "Norm is:  1.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1086ded50>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHRBJREFUeJzt3X2UXHWd5/H3hw4ZIREiJCSQNIRAwER50mNQQaZUVno8\naHR2VwjquI4zm9mduHt21EV0jjQy7sg5uKOzKMaZjMfj8RiddYAMCuFBanyAReKEhIE0JEJId+cR\nkIfwmE6++8etxErR6arqrqr7UJ/XOXXoW/fWvV/bzqd//b33/q4iAjMzK4bD0i7AzMxax6FuZlYg\nDnUzswJxqJuZFYhD3cysQBzqZmYFUjfUJfVJGpC0UdLlo6x/naQbJK2TdK+kN7SnVDMzq2fMUJfU\nA1wH9AELgSWSFtRs9jngXyPiLOCPgK+1o1AzM6uv3kh9EbApIjZHxB5gJbC4ZpsFwF0AEfEwMFfS\njJZXamZmddUL9dnAYNXyUOW9auuAPwSQtAg4CZjTqgLNzKxx9UK9kTkEvgxMk7QWWAasBfZOtDAz\nM2vepDrrh4HequVektH6ARHxHPDH+5clPQY8WrsjSZ5kxsxsHCJCjW5bb6S+Bpgvaa6kycAlwKrq\nDSQdXVmHpD8F/iUidh+iML9a9LryyitTryHvr1deCRYtCl7zmiv5whfSr6coL/9stvbVrDFDPSJG\nSFoqq4GHgB9ExAZJSyUtrWy2EHhA0gBwEfDfm67CLAXf/CZMmwZ//ufwjW/Axo1pV2Q2cfXaL0TE\nLcAtNe8tr/r6HuD01pdm1j4RcP318K1vwR13wMc/DitWwJe/nHZlZhNTN9Qtm0qlUtol5Nq6dfDy\ny3DeeTAyUuK44+Cii+Cv/xrUcPfSRuOfzXRpPD2bcR1Iik4dy6yea6+Fxx6Dr389WY6AefPg5pvh\nDb4n2jJEEtHCE6VmhXTXXfCud/1uWYL3vAduuy29msxawaFuXWfPHvjFL6C2S3DBBXD33amUZNYy\nDnXrOuvWwUknwbHHHvz+W94C992XTk1mreJQt66zdi286U2vfv/UU+Hpp2HXrs7XZNYqDnXrOvff\nD2ef/er3DzsM3vxmWLOm8zWZtYpD3brO/ffDOeeMvs4tGMs7h7p1lX37YP16OOus0defcQY8+GBn\nazJrJYe6dZVHH01OkE6bNvr6BQtgw4bO1mTWSg516yoDA0lwH8rppydzwIyMdK4ms1ZyqFtX2bgR\n5s8/9PopU2DWrORuU7M8cqhbV3nkETjttLG3cQvG8syhbl2l3kgdHOqWbw516yqNjNRPO81zq1t+\nOdSta7z4YnK36Iknjr3dvHnuqVt+1Q11SX2SBiRtlHT5KOunS7pV0v2S/k3Sf2pLpWYTtGkTzJ0L\nPT1jb3fyycmlj2Z5NGaoS+oBrgP6SB5bt0RS7QVhy4C1EXE2UAK+IskP37DM2bixfusFkpH81q3J\nbI5meVNvpL4I2BQRmyNiD7ASWFyzzTbgqMrXRwFPVp5tapYpjzxS/yQpwOTJyWWNg4Ptr8ms1eqF\n+myg+kd7qPJetb8D3iBpK7AOP3jaMmrz5qS10gj31S2v6rVJGnn+3OeA+yOiJOkU4HZJZ0XEc7Ub\n9vf3H/i6VCr5WYbWUY8/Dhdf3Ni2+/vq7353e2syq1UulymXy+P+fL1QHwZ6q5Z7SUbr1d4OfAkg\nIn4j6THgdOBVE5hWh7pZp23ZkjwcoxEeqVtaage8V111VVOfr9d+WQPMlzRX0mTgEmBVzTYDwIUA\nkmaSBLqvHbBMiUhG6vUuZ9zPV8BYXo05Uo+IEUnLgNVAD7AiIjZIWlpZvxz4X8C3Ja0j+SXxPyPi\nqTbXbdaUp56CSZPg6KMb294jdcsrRTTSNm/BgaTo1LHMaq1dCx/7WDKXeiO2bUuejrRjR3vrMqtH\nEhGhRrf3HaXWFR5/vPF+OsDMmcnzSl9+uX01mbWDQ926QjMnSSF5Xunxx8PwcPtqMmsHh7p1hWZO\nku43Zw4M1V7rZZZxDnXrCs22XwB6e31XqeWPQ926wpYtHqlbd3CoW1cYHExG3s3wSN3yyKFuhTcy\nAk88kUzS1QyP1C2PHOpWeDt2wPTpyc1HzfBI3fLIoW6FNzwMJ5zQ/Oc8Urc8cqhb4W3dCrNrJ4xu\ngG9AsjxyqFvhjXek7huQLI8c6lZ44x2pg1swlj8OdSu8rVvHN1KHZKS+bVtr6zFrJ4e6Fd7w8PhH\n6g51yxuHuhWeR+rWTRzqVnjjPVEKSahv3draeszaqW6oS+qTNCBpo6TLR1n/aUlrK68HJI1Imtae\ncs2a88IL8OKLcOyx4/u8R+qWN2OGuqQe4DqgD1gILJG0oHqbiLg2Is6JiHOAK4ByRDzdroLNmrFt\nWxLMavi5MQdzqFve1BupLwI2RcTmiNgDrAQWj7H9ZcD3W1Wc2URN5CQpONQtf+qF+mygevaLocp7\nryLpSOAi4EetKc1s4iZykhSSts3u3fDSS62ryayd6k1x1MyTot8H/GKs1kt/f/+Br0ulEqVSqYnd\nmzVvoiP1ww5LZnfcvh3mzm1ZWWaHVC6XKZfL4/58vVAfBqpnoe4lGa2P5lLqtF6qQ92sE/b31Cdi\nfwvGoW6dUDvgveqqq5r6fL32yxpgvqS5kiYDlwCrajeSdDRwAXBTU0c3a7Pt25ufR72W++qWJ2OO\n1CNiRNIyYDXQA6yIiA2SllbWL69s+gFgdUS82NZqzZq0Y0cy2+JEONQtT+o+NiAibgFuqXlvec3y\nd4DvtLY0s4lrVaj7BiTLC99RaoXmkbp1G4e6FdbICDz1FMyYMbH9nHCCQ93yw6FuhfXEE/C61zX/\nbNJaHqlbnjjUrbBa0XqB5OoZh7rlhUPdCqtVoT5jBjz5JOzdO/F9mbWbQ90KqxXXqAMcfjhMm5a0\nc8yyzqFuhdWqkTok+9mxozX7Mmsnh7oVlkPdupFD3QqrVe0XSEJ9587W7MusnRzqVlgeqVs3cqhb\nYTnUrRs51K2wHOrWjRzqVkitmiJgP4e65YVD3QqpVVME7OdQt7xwqFshtbL1Ag51yw+HuhVSKy9n\nBDjuONi1C/bta90+zdqhbqhL6pM0IGmjpMsPsU1J0lpJ/yap3PIqzZrU6pH67/0eTJkCv/1t6/Zp\n1g5jdhwl9QDXAReSPIT6PkmrImJD1TbTgK8DF0XEkKTp7SzYrBGtDnX4XQvm2GNbu1+zVqo3Ul8E\nbIqIzRGxB1gJLK7Z5jLgRxExBBARnvbIUtfOUDfLsnqhPhsYrFoeqrxXbT5wjKS7JK2R9NFWFmg2\nHq3uqYND3fKh3gVf0cA+DgfeBLwbOBK4R9L/i4iNtRv29/cf+LpUKlEqlRou1KwZHqlbXpXLZcrl\n8rg/Xy/Uh4HequVektF6tUHgiYh4EXhR0s+As4AxQ92snRzqlle1A96rrrqqqc/Xa7+sAeZLmitp\nMnAJsKpmm5uA8yX1SDoSOBd4qKkqzFpsxw63X6w7jTlSj4gRScuA1UAPsCIiNkhaWlm/PCIGJN0K\nrAf2AX8XEQ51S83+KQKmt/g6LIe65YEiGmmbt+BAUnTqWNbdtm+HM89s/fzn994Ly5bBffe1dr9m\nY5FERKjR7X1HqRVOO/rp4JG65YND3QqnHZczwu+efuQ/OC3LHOpWOO0aqR9xBEyeDM880/p9m7WK\nQ90Kp12hDm7BWPY51K1w2nE5434Odcs6h7oVzvbtHqlb93KoW+G4/WLdzKFuheNQt27mULfCcahb\nN3OoW6Hs3ZtMETBjRnv271C3rHOoW6E88QRMmwaT6s0/Ok4Odcs6h7oVSjtbL+BQt+xzqFuhdCrU\nPVWAZZVD3Qql3aE+dSocdhjs3t2+Y5hNhEPdCqXdoQ7J/rdvb+8xzMarbqhL6pM0IGmjpMtHWV+S\n9IyktZXXX7anVLP6OhXq7qtbVo15jYCkHuA64EKS55XeJ2lVRGyo2fRfIuL9barRrGE7dsCCBe09\nhkPdsqzeSH0RsCkiNkfEHmAlsHiU7Rp+KodZO3mkbt2uXqjPBgarlocq71UL4O2S1kn6iaSFrSzQ\nrBmdCPVZs9xTt+yqF+qNXLj1r0BvRJwF/B/gxglXZTZOHqlbt6t3390w0Fu13EsyWj8gIp6r+voW\nSd+QdExEPFW7s/7+/gNfl0olSqXSOEo2G92+fbBrFxx3XHuP41C3diqXy5TL5XF/XjHGXRSSJgEP\nA+8GtgK/ApZUnyiVNBPYGREhaRHww4iYO8q+YqxjmU3UE0/Aaaclc7+00y9/CZ/+NNxzT3uPYwYg\niYho+LzlmCP1iBiRtAxYDfQAKyJig6SllfXLgf8A/BdJI8ALwKXjrt5sAjrReoGkp+6RumXVmCP1\nlh7II3Vrs5/+FL74RZjAX64N2b07afE8/zzI131ZmzU7UvcdpVYYnRqpT52ahLmnCrAscqhbYXQq\n1MEnSy27HOpWGJ0MdffVLasc6lYYnR6p+wYkyyKHuhXGzp3tv0Z9P7dfLKsc6lYY7qmbOdStQBzq\nZg51K4iIpP3SyROl7qlbFjnUrRCeeQYmT4YjjujM8TxSt6xyqFshdLL1Ag51yy6HuhVCWqHumS8s\naxzqVgidDvWpU5P/eqoAyxqHuhVCp0NdcgvGssmhboXQ6VAHh7plk0PdCsGhbpZwqFshpBHqntTL\nsqhuqEvqkzQgaaOky8fY7i2SRiT9YWtLNKsvrZG6b0CyrBkz1CX1ANcBfcBCYImkBYfY7hrgVsDP\ngrGOc/vFLFFvpL4I2BQRmyNiD7ASWDzKdp8E/i+wq8X1mdUVkYyYZ83q7HEd6pZF9UJ9NjBYtTxU\nee8ASbNJgv76ylu+HcM66pln4PDDYcqUzh7XoW5ZVC/UGwnorwKfrTxVWrj9Yh22bRscf3znj+tJ\nvSyLJtVZPwz0Vi33kozWq70ZWKnkserTgT+QtCciVtXurL+//8DXpVKJUqnUfMVmNdIKdY/UrR3K\n5TLlcnncn1eMMXmFpEnAw8C7ga3Ar4AlEbHhENt/G/jniPinUdbFWMcyG6/vfQ9uvhm+//3OHjci\nafns3Pm7aQPMWk0SEdFwB2TM9ktEjADLgNXAQ8APImKDpKWSlk6sVLPWSGuk7qkCLIvqtV+IiFuA\nW2reW36IbT/eorrMGpZWqMPv+uqnnJLO8c1q+Y5Sy700Q/3445Pjm2WFQ91yL81QP+EEh7pli0Pd\nci/tkfrWrekc22w0DnXLvbRH6g51yxKHuuXaCy/AK6/AtGnpHN+hblnjULdc27YtuQJFKd3H7FC3\nrHGoW66l2XoBnyi17HGoW66lHerHHAPPPw8vvpheDWbVHOqWa9u3pxvqkq9Vt2xxqFuupT1SB/fV\nLVsc6pZrDnWzgznULdcc6mYHc6hbrmUh1N1TtyxxqFuubd2afqh7pG5Z4lC33HrpJXj2WZgxI906\nHOqWJQ51y62tW5NAPSzln2KHumVJ3X8OkvokDUjaKOnyUdYvlrRO0lpJv5b0rvaUanawoSGYPTvt\nKhzqli1jPvlIUg9wHXAhyUOo75O0quYZpXdExE2V7c8AbgBObVO9ZgcMDcGcOWlXkUwm9soryZ2l\nU6akXY11u3oj9UXApojYHBF7gJXA4uoNIuL5qsWpwBOtLdFsdMPD2Qh131VqWVIv1GcDg1XLQ5X3\nDiLpA5I2kDzL9L+1rjyzQ8vKSB2SFszwcNpVmNV/8HQ0spOIuBG4UdI7gO8Cp4+2XX9//4GvS6US\npVKpoSLNRjM0BOefn3YViTlzknrMJqpcLlMul8f9eUUcOrclvRXoj4i+yvIVwL6IuGaMz/wGWBQR\nT9a8H2Mdy6xZ554LX/0qvO1taVcCn/kMHHssfPazaVdiRSOJiGj4iQH12i9rgPmS5kqaDFwCrKo5\n4ClS8ogCSW8CqA10s3bIUvultxcGB+tvZ9ZuY7ZfImJE0jJgNdADrIiIDZKWVtYvB/498EeS9gC7\ngUvbXLMZIyOwa1fy1KMs6O2FO+5IuwqzOu2Xlh7I7RdroaEhWLQoO9eH//rX8Cd/AmvXpl2JFU2r\n2y9mmZSl1gskI/UtW9KuwsyhbjmVtVCfMSO5+eiFF9KuxLqdQ91yKStTBOwnJb9kfLLU0uZQt1zK\nWqiDr4CxbHCoWy5t2QInnZR2FQdzqFsWONQtlzZvhrlz067iYCee6FC39DnULZeyGOoeqVsWONQt\nd154IXni0cyZaVdyMIe6ZYFD3XLn8ceTVkfaTzyq5VC3LMjYPwuz+h5/PHsnSSH5RbNlC/jGaUuT\nQ91yJ4v9dICjj4bDD4cnPZ2dpcihbrnz+OPZDHWAk0+GRx9NuwrrZg51y53Nm7PZfgGYNw8eeyzt\nKqybOdQtdzxSNzs0h7rljkfqZofmULdceeml5ETkCSekXcno5s3zSN3S1VCoS+qTNCBpo6TLR1n/\nYUnrJK2X9EtJZ7a+VLPkOvA5c6CnJ+1KRnfyyR6pW7rqhrqkHuA6oA9YCCyRtKBms0eBCyLiTOBq\n4FutLtQMksDMaj8dkrbQ0FDyuD2zNDQyUl8EbIqIzRGxB1gJLK7eICLuiYhnKov3Ahl6fIEVyaZN\ncOqpaVdxaJMnJ9MXDA2lXYl1q0ZCfTZQffPzUOW9Q/kE8JOJFGV2KFkPdXBf3dI1qYFtGr7pWdI7\ngT8GzhttfX9//4GvS6USpVKp0V2bAUmoX3BB2lWMzX11m4hyuUy5XB735xV1JqqQ9FagPyL6KstX\nAPsi4pqa7c4E/gnoi4hNo+wn6h3LrJ4FC+Af/xHe+Ma0Kzm0q69OrtL50pfSrsSKQBIRoUa3b6T9\nsgaYL2mupMnAJcCqmoOeSBLoHxkt0M1aYe/eZAQ8b17alYztlFOSvyjM0lC3/RIRI5KWAauBHmBF\nRGyQtLSyfjnwBeB1wPWSAPZExKL2lW3daGgIpk+HI49Mu5KxnXYaPPxw2lVYt6rbfmnZgdx+sQm6\n8074q7+Cu+5Ku5KxPfssHH88PPdc9uZ8t/xpR/vFLBM2bsz+lS8ARx2VTMM7PJx2JdaNHOqWG3kJ\ndYDTT3cLxtLhULfcGBiA178+7Soa4766pcWhbrnx0EOwcGHaVTTGI3VLi0PdcuH552H79uTGnjxw\nqFtaHOqWCwMDSUtjUiP3QGeAQ93S4lC3XNiwIT+tF0j+oti5M7ms0ayTHOqWC3nqp0My3/uCBfDg\ng2lXYt3GoW65kLdQBzjjDHjggbSrsG7jULdccKibNcahbpm3e3cy70tebjza78wzHerWeQ51y7z1\n65NR+uGHp11Jc844I6ndUx5ZJznULfPWroVzzkm7iubNnJlM6LVtW9qVWDdxqFvm5TXUpaQFs359\n2pVYN3GoW+blNdQB3vxmWLMm7SqsmzjULdP27EluPDrzzLQrGZ9Fi+BXv0q7CusmDYW6pD5JA5I2\nSrp8lPWvl3SPpJckfar1ZVq3eughOOkkmDIl7UrG5y1vgfvu88lS65y6oS6pB7gO6AMWAkskLajZ\n7Engk8C1La/Qutq99ybBmFcnnpg8W9UPzLBOaWSkvgjYFBGbI2IPsBJYXL1BROyKiDXAnjbUaF3s\n7rvhvPPSrmL8pOSXklsw1imNhPpsYLBqeajynlnb3X03vP3taVcxMe6rWyc1MpFpy7qB/f39B74u\nlUqUSqVW7doKaOfO5PWGN6RdycS87W1w9dVpV2F5US6XKZfL4/68os4ZHElvBfojoq+yfAWwLyKu\nGWXbK4HdEfGVUdZFvWOZVbvpJrj+erj11rQrmZjdu2HWrOQX1JFHpl2N5Y0kIkKNbt9I+2UNMF/S\nXEmTgUuAVYc6fqMHNqvnl7/Mf+sFYOpUOOssuOeetCuxblA31CNiBFgGrAYeAn4QERskLZW0FEDS\nLEmDwP8A/lLSFklT21m4Fd8dd8C73pV2Fa1RKsEE/qI2a1jd9kvLDuT2izVh587k8XW7duVvIq/R\n3H47fPGL8POfp12J5U072i9mHXf77fDOdxYj0CFpI61bB88+m3YlVnQOdcuk1avhoovSrqJ1pkyB\n889P/neZtZND3TJn71647TZ4z3vSrqS13v9+WHWoSwzMWsShbpnzi18klwDOm5d2Ja118cXwk5/A\nyEjalViROdQtc37wA7jkkrSraL05c+Dkk32y1NrLoW6ZMjICP/oRfOhDaVfSHpdeCt/9btpVWJE5\n1C1TfvpT6O2FU05Ju5L2+MhH4IYbkrtMzdrBoW6Zsnw5fOITaVfRPrNmwTvekfw1YtYOvvnIMmN4\nGN74RtiyBV772rSraZ+bboIvfSmZK16eWMPq8M1Hllt/8zfw0Y8WO9AB3vc+eO65pNVk1moeqVsm\n7NoFp58O69cnV4kU3Xe+k7wc7FaPR+qWS1deCR/+cHcEOsBll8HQEPz4x2lXYkXjkbql7te/hve+\nFzZsgGOOSbuazrntNvizP4MHH4Qjjki7Gssqj9QtV557DpYsga99rbsCHZJpEM49Fz71qbQrsSLx\nSN1Ss3dvcjPOUUfBihVpV5OOZ55JHkx9xRXw8Y+nXY1lUctH6pL6JA1I2ijp8kNs87eV9eskndNM\nwdad9uxJrkd/+mn4+tfTriY9Rx8NN94In/uc7zS11hgz1CX1ANcBfcBCYImkBTXbvBc4NSLmA/8Z\nuL5NtVqViTyYNm2Dg3DhhcmDMG68EV7zmnTrSft7uXAh3HlnEuyf+Qy88kqq5UxY2t/PbldvpL4I\n2BQRmyNiD7ASWFyzzfuB7wBExL3ANEkzW16pHSSP/3AGB+Hzn4ezz076yTffnMwznrYsfC8XLoS1\na2HTJliwAL79bXj++bSrGp8sfD+72aQ662cDg1XLQ8C5DWwzB9gx4eost156CTZvhkceSa5uufPO\n5OqWyy6DNWuS2QrtYNOnJ/PC3HUXXHst/MVfJE9/uuCC5MHVp5wCs2dDT0/alVqW1Qv1Rs9s1jbx\nR/3cxRdXbVCzRb1lf+bg9wYHk6foZKm2l1+G3/42ee3bByedlDxn9IwzkuvQzz/fl+414p3vTF7b\ntyc3J/3sZ0nY/+Y3yXtTp8K0ackJ5smTYdKkV7/Gmn7gUOvG85nR1j38cPKL3NIx5tUvkt4K9EdE\nX2X5CmBfRFxTtc03gXJErKwsDwC/HxE7avblS1/MzMahmatf6o3U1wDzJc0FtgKXAEtqtlkFLANW\nVn4JPF0b6M0WZWZm4zNmqEfEiKRlwGqgB1gRERskLa2sXx4RP5H0XkmbgOcBX21rZpaSjt18ZGZm\n7dfWaQIk/UdJD0raK+lNNeuuqNywNCCpYM+Nbz9J/ZKGJK2tvPrSrimPGrm5zhonabOk9ZWfyV+l\nXU+eSPoHSTskPVD13jGSbpf0iKTbJE2rt592z/3yAPBB4GfVb0paSNKfX0hyY9M3JHkemuYE8L8j\n4pzK69a0C8qbRm6us6YFUKr8TC5Ku5ic+TbJz2K1zwK3R8RpwJ2V5TG1NUgjYiAiHhll1WLg+xGx\nJyI2A5tIbnSy5vjk88Q0cnOdNc8/l+MQET8Hflvz9oGbOyv//UC9/aQ1Oj6B5Cal/YZIbmKy5nyy\nMt/Oikb+LLNXGe3GOf8cTkwAd0haI+lP0y6mAGZWXU24A6h7t369SxrrknQ7MGuUVZ+LiH9uYlc+\nY1tjjO/t50nm2PliZflq4CtAgR/Z3Bb+mWu98yJim6QZwO2SBiojUJugiIhG7veZcKhHxL8bx8eG\ngd6q5TmV96xKo99bSX8PNPML1BK1P4e9HPwXpDUpIrZV/rtL0g0kLS6H+vjtkDQrIrZLOh7YWe8D\nnWy/VPfZVgGXSpos6WRgPuAz5U2o/B+83wdJTkpbcw7cXCdpMsnJ+1Up15Rbko6U9NrK11OA9+Cf\ny4laBXys8vXHgBvrfWDCI/WxSPog8LfAdODHktZGxB9ExEOSfgg8BIwA/9VP0GjaNZLOJmkhPAYs\nTbme3DnUzXUpl5VnM4EblEwGMwn4XkTclm5J+SHp+8DvA9MlDQJfAL4M/FDSJ4DNwIfq7sdZamZW\nHL423MysQBzqZmYF4lA3MysQh7qZWYE41M3MCsShbmZWIA51M7MCcaibmRXI/wcNLE0t6kB/jQAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1085d0c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from math import pi\n",
    "root_pi = np.sqrt(pi)\n",
    "# Define a gaussian with proper normalisation as our test function\n",
    "def gauss(x,alpha):\n",
    "    return np.sqrt(alpha / (root_pi))* np.exp(-0.5 * alpha**2 * x**2)\n",
    "# We can also write the second derivative function easily enough\n",
    "def d2gauss(x,alpha):\n",
    "    return np.sqrt(alpha / (root_pi))* alpha*alpha * (alpha*alpha*x*x - 1) * np.exp(-0.5 * alpha**2 * x**2)\n",
    "\n",
    "# Declare space for the potential and call the routine\n",
    "omega = 2.0 # Set the frequency\n",
    "VQHO = np.linspace(0.0,width,num_x_points)\n",
    "square_well_potential(x,VQHO,omega*omega)\n",
    "print \"From theory, E min, alpha: \",0.5*omega,np.sqrt(omega)\n",
    "psi = gauss(x,np.sqrt(omega))\n",
    "# Check that I've normalised the gaussian correctly\n",
    "print \"Norm is: \",integrate_functions(psi,psi,num_x_points,dx)\n",
    "pl.plot(x,psi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the trial function, with the parameter $\\alpha$ as our variational parameter, we can vary the energy and look for the minimum.  We'll do this in a **very** crude way, just scanning values of $\\alpha$.  There are much better approaches, though, which we will address in a later notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy at optimum alpha is:  1.0\n",
      "Min E and alph:  1.0000178072 1.41\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10868d210>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAERCAYAAACdPxtnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XncnfOd//HXWyKbLVREYmkstdNYR0O5FR1FlWotY4lp\nR2vasXWYYtrpPdqJbRr16wxGLA1F1ZYHVSUiZ4YqGg0iobFXShKSEBGy3Z/fH99zc7vdue9z3znX\nfZ3rnPfz8TiP+yzXuc7nWK7P+XxXRQRmZtaYVss7ADMzy4+TgJlZA3MSMDNrYE4CZmYNzEnAzKyB\nOQmYmTWwzJOApD6Spkq6u/y4WdKs8nNTJR2UdQxmZtaxvr3wGacDM4C1yo8DGBsRY3vhs83MrBOZ\nVgKSNgYOBq4G1Pp0m/tmZpajrJuDLgXOBlraPBfAqZKeknSNpMEZx2BmZiuRWRKQdCgwNyKm8vFf\n/lcAmwEjgTeAn2YVg5mZdU5ZrR0kaQxwArAcGACsDdweESe2OWYEcHdE7NjB+72okZlZD0RExU3u\nmVUCEXFeRGwSEZsBxwAPRsSJkoa1OewIYFon56jb249+9KPcY/B38/fz96u/W3f1xuggSM1BrdFd\nLOmz5ccvA9/upRjMzKydXkkCEVECSuX7J/TGZ5qZWdc8YzgnTU1NeYeQmXr+buDvV3T1/v26K7OO\n4VUlKWo1NjOzWiWJqIWOYTMzq31OAmZmDcxJwMysgdV8Enj8cWhp6fo4MzPrvprvGF53XXjhBfjU\np/KOyMys9tVdx/CGG8KcOXlHYWZWn2o+CQwdCrNn5x2FmVl9qvkkMGQIvPlm3lGYmdWnmk8Ca6wB\n77+fdxRmZvWp5pPAoEGweHHeUZiZ1aeaTwIDBzoJmJllpeaTgCsBM7PsOAmYmTUwJwEzswaWeRKQ\n1EfSVEl3lx+vJ2mipJmS7pc0uLP3OwmYmWWnNyqB04EZfLS95DnAxIjYCphUfrxSgwZ5iKiZWVYy\nTQKSNgYOBq4m7TMMcBgwvnx/PHB4Z+dwJWBmlp2sK4FLgbOBtuuADo2I1tWA5gBDOzuBk4CZWXYy\n22he0qHA3IiYKqmpo2MiIiStdBnT5uZmXn0Vnn0WSqUm7w1qZtZOqVSiVCr1+P2ZLSUtaQxwArAc\nGACsDdwB7A40RcRsScOAyRGxTQfvj4hgyhT49rfhiScyCdPMrK7UzFLSEXFeRGwSEZsBxwAPRsQJ\nwF3A6PJho4EJnZ3HzUFmZtnpzXkCrSXHhcCBkmYCXyg/XiknATOz7NT8zmKzZ8PIkd5TwMysEjXT\nHFQt/frBkiV5R2FmVp9qPgn07w9Ll+YdhZlZfar5JOBKwMwsOzWfBPr2hZYWWLEi70jMzOpPzScB\nKVUDbhIyM6u+mk8CkPoF3CRkZlZ9hUgCrgTMzLJRiCTgSsDMLBuFSAKuBMzMslGIJOBKwMwsG4VJ\nAq4EzMyqrxBJwBPGzMyyUYgk4ErAzCwbhUgCrgTMzLJRiCTgjmEzs2wUIgl4iKiZWTYyTQKSBkh6\nTNKTkmZIuqD8fLOkWZKmlm8HdXYeVwJmZtnom+XJI+IDSftFxGJJfYGHJe1N2mpybESMreQ8rgTM\nzLKReXNQRLTuENwP6AMsKD+uePszVwJmZtnIPAlIWk3Sk8AcYHJETC+/dKqkpyRdI2lwZ+fwEFEz\ns2xk2hwEEBEtwEhJ6wD3SWoCrgDOLx/yY+CnwDfbv7e5uRmAJ56AiCagKetwzcwKpVQqUSqVevx+\nRUT1ounqw6QfAu9HxH+2eW4EcHdE7Nju2GiN7ZxzYPDg9NfMzFZOEhFRcXN71qOD1m9t6pE0EDgQ\nmCppwzaHHQFM6+w8nixmZpaNrJuDhgHjJa1GSjg3RMQkSddLGkkaJfQy8O3OTtK/Pyxe3NkRZmbW\nE1kPEZ0G7NLB8yd25zz9+sGCBV0fZ2Zm3VOIGcMeImpmlo1CJAFPFjMzy0YhkoArATOzbBQiCbgS\nMDPLRiGSgCsBM7NsFCYJuBIwM6u+QiQBTxYzM8tGIZKAKwEzs2wUIgm4EjAzy0YhkoArATOzbBQi\nCbgSMDPLRiGSgCsBM7NsFCIJuBIwM8tGIZKAKwEzs2wUIgm4EjAzy0YhkoArATOzbGSWBCQNkPSY\npCclzZB0Qfn59SRNlDRT0v2t2092xpWAmVk2Mt1oXtKgiFgsqS/wMHAWcBjwVkRcLOn7wLoR8Ykt\n5NtuNB8Bq60GK1akv2Zm1rGa2mg+Ilp3Bu4H9AEWkJLA+PLz44HDuzqP5OWkzcyykGkSkLSapCeB\nOcDkiJgODI2IOeVD5gBDKzmXl5M2M6u+rDeabwFGSloHuE/Sfu1eD0krbY9qbm7+8L7UxNKlTRlF\namZWTKVSiVKp1OP3Z9on8LEPkn4IvA/8A9AUEbMlDSNVCNt0cHy0jW2jjeCxx2DjjXslXDOzQqqZ\nPgFJ67eO/JE0EDgQmArcBYwuHzYamFDJ+dwnYGZWfVk2Bw0DxktajZRsboiISZKmAr+W9E3gFeCo\nSk7mPgEzs+rLLAlExDRglw6enw8c0N3zuRIwM6u+woy6dyVgZlZ9hUkCrgTMzKqvMEnAlYCZWfUV\nJgm4EjAzq77CJAFXAmZm1VeYJOBKwMys+gqTBFwJmJlVX2GSgCsBM7PqK0wScCVgZlZ9hUoCrgTM\nzKqrMEnAW0yamVVfYZKAKwEzs+orTBJwJWBmVn2FSQKuBMzMqq8wScCVgJlZ9WW90fwmkiZLmi7p\nGUmnlZ9vljRL0tTy7aCuzuVKwMys+jLdaB5YBpwZEU9KWhN4QtJEIICxETG20hO5EjCzapoyBXbb\nLe8o8tdlJSBprKTte3LyiJgdEU+W7y8CngU2aj11d87lSsDMqiECxoyBr38d5s/PO5r8VdIc9Cxw\nlaTHJZ0iaZ2efJCkEcDOwKPlp06V9JSka1o3pO+MKwEzW1UrVsB3vgO33gqPPALrrZd3RPnrMglE\nxLiI2As4ERgBTJN0k6T9Kv2QclPQbcDp5YrgCmAzYCTwBvDTrs7hSsDMVsXixXDkkfDCC/C//wvD\nhuUdUW2oqE9AUh9gG2Bb4E3gKeB7kk6JiKO7eO/qwO3ALyNiAkBEzG3z+tXA3R29t7m5+cP7Awc2\nsWRJUyXhmpl9zFtvwZe/DFtuCb/+dWpZqBelUolSqdTj9ysiOj9AuhT4MvAgcHVEPN7mtT9HxNad\nvFfAeGBeRJzZ5vlhEfFG+f6ZwO4R8Xft3httY3vwQfjJT9JfM7NKvfQSfOlL8NWvpr4Adas3sngk\nEREVf8tKKoGngR9ExHsdvPY3Xbx3L+B44GlJU8vPnQccK2kkaZTQy8C3uwrCfQJm1l2PPgpHHAE/\n+AF897t5R1ObKk0CW+vj6fMd4NWIeLuzN0bEw3Tc73BvxRGWuU/AzLrj9tvhlFPguuvg0EPzjqZ2\nVZIE/hvYlZQMAHYEpgPrSPrHiLgvq+DaciVgZpWIgLFj4dJL4b77YJdd8o6otlUyRPR1YGRE7BoR\nu5JG9LwEHAhcnGVwbXlTGTPryvLlqdnnF79IQ0CdALpWSSWwdURMb30QETMkbRMRL0rqvFe5igYM\ngA8+6K1PM7OiWbQIjj4ali2Dhx+GdXo0o6nxVFIJTJd0haR9JTVJuhyYIak/aVmIXjFwoJOAmXXs\n9ddhn33S2P977nEC6I5KksBo4EXgDOB0UlPQaFIC+EJ2oX2cKwEz68i0afC5z8HXvgbjxsHqq+cd\nUbF0Ok9AUl9gYkRUPDu4WtrPE1iyBNZayyOEzOwj998Pxx8Pl10Gxx6bdzS1obvzBDqtBCJiOdBS\nydo+WevXL3X6rFiRdyRmVgv++7/hxBPTUFAngJ6rpGP4PdJ6QRPL9wEiIk7LLqxPklKT0JIlMGhQ\nb36ymdWSZcvgjDOgVEojgDbfPO+Iiq2SJHBH+dbaNqM293tVa7+Ak4BZY1qwIC0B3a9fSgDuAF51\nXSaBiPiFpEHAphHxXC/EtFIDBsD77+cZgZnlZebMtAjcwQfDf/4n9OmTd0T1oZJNZQ4DpgK/Kz/e\nWdJdWQfWEY8QMmtMkybB5z8PZ52VZgI7AVRPJUNEm0kLxS0AiIipQC6tcJ4rYNZ4rrgCjjsObrkF\nTj4572jqTyV9Assi4u12C8i1ZBRPp1wJmDWO5cvhzDPhgQfSDOAtt8w7ovpUSRKYLuk4oK+kzwCn\nAY9kG1bH3Cdg1hgWLIBjjkmjAv/wBxic+yD1+lVJc9CpwPbAEuBmYCFp9nCvcyVgVv+mT4c99oBt\nt4Xf/MYJIGuVjA56j7QRzHnZh9M59wmY1bcJE1K7/yWXwEkn5R1NY+gyCUjaGjiLtMl86/EREV2u\nGyRpE+B6YAPS3IKrIuL/SVoPuAX4NPAKcFRXG9SAm4PM6lVLC/z7v6cNYH77W9h997wjahyV9Anc\nClwBXA20LtpQ6WSxZcCZEfGkpDWBJ8ozj/+etCbRxZK+D5xTvnXKzUFm9WfhQjjhBJg/H/74Rxg6\nNO+IGkulo4Ou6MnJI2I2MLt8f5GkZ4GNgMOAfcuHjQdKVJAE3BxkVl9mzoSvfAX22w9uvTXNBLbe\nVUnH8N2SvitpmKT1Wm/d/SBJI4CdgceAoRExp/zSHKCi3O9KwKx+3HMP7L03fO97cPnlTgB5qaQS\nOInU/HNWu+c3q/RDyk1BtwOnR8S7beccRERUukOZ+wTMii8CLrggrQI6YQKMGpV3RI2tktFBI1bl\nAyStTkoAN0TEhPLTcyRtGBGzJQ0D5nb03ubm5g/vNzU1MWBAkysBswJbuBC+8Q2YNSu1/w8fnndE\nxVcqlSiVSj1+/0o3lZH0LxFxcfn+1yPi1javjYmILoeMKv3kHw/Mi4gz2zx/cfm5iySdAwyOiHPa\nvTfaxzZmTNpHdMyYyr+gmdWG6dPhyCOhqQl+9rNU2Vv1VXNTmbbbNLS/4H+pwvPvBRwP7Cdpavl2\nEHAhcKCkmaQtKi+s5GTuEzArpptvThf/c8+FK690AqgllfQJ9FhEPMzKE80B3T2f+wTMimXpUvjn\nf4Z7701rAH32s3lHZO1lmgSqzZWAWXHMmpU2gNlgA5gyxcs/1KrOmoN2kvSupHeBHVvvtz7upfg+\nxvMEzIph0qQ06/fww+HOO50AatlKK4GIqLltGwYMgMWL847CzFampQUuvBB+/nO48Ub4QpeLy1je\nCtUctMYa8N57XR9nZr1vwQIYPRreeis1/2y0Ud4RWSUqmTFcM5wEzGrTo4/CzjvD5ptDqeQEUCSu\nBMysx1paYOxYuPhiuOqq1AdgxeIkYGY9Mm9eav6ZNy/N/v30p/OOyHrCzUFm1m0PP5yaf7bbDv7v\n/5wAiqxQlcCaazoJmOWppQUuugguuwyuuQYOOSTviGxVFSoJrLFGGiIakTagNrPeM3du2vxl8eI0\n+mfjjfOOyKqhUM1BffrA6qt7wphZb3vwQdhlF9htN5g82QmgnhSqEoCP+gUGDsw7ErP6t3Qp/PCH\n8Mtfpv1/v/jFvCOyaitsElh//bwjMatvM2fC3/0dDBsGTz4JQ4bkHZFloVDNQeARQmZZi4Crr4a9\n9kobwNx1lxNAPStsJWBm1Td/PnzrW/D882nm7/bb5x2RZa2QlcCiRXlHYVZ/Jk9O6/1vsgk89pgT\nQKPINAlIulbSHEnT2jzXLGlWu53GKuZKwKy6li5NO34ddxyMGweXXuqdvxpJ1s1B1wE/B65v81wA\nYyNibE9O6CRgVj0zZsCJJ8LQoanzd4MN8o7IelumlUBEPAQs6OClHk/1chIwW3UtLekX/z77wMkn\nw29+4wTQqPLqGD5V0onAFOCfI+LtSt/oJGC2al55BU46CZYvT23/W2yRd0SWpzySwBXA+eX7PwZ+\nCnyzowObm5s/vN/U1ERTU5OTgFkPRcC118I558DZZ6cN4PvU3P6B1l2lUolSqdTj9ysiqhdNRx8g\njQDujohP7EvcxWvRUWznn586sn7yk6qHala3Zs9OzT6zZsH118OOuewSbr1BEhFRcZN7rw8RlTSs\nzcMjgGkrO7Yja68NCxdWNyazenbbbTByZBr++dhjTgD2cZk2B0m6GdgXWF/Sa8CPgCZJI0mjhF4G\nvt2dc66zDrzzTtVDNas78+bBaaelDV8mTIA998w7IqtFmSaBiDi2g6evXZVzDh7sJGDWldtvh1NP\nhaOOgqlT04AKs44UbtkIVwJmKzdnDvzTP8G0aakZaNSovCOyWle4ZSOcBMw+KQJuvBF22gm23DJN\n/HICsEq4EjAruL/+FU45JY3/v+eetPGLWaUKWQm8XfHUMrP6FZH2+R05EnbdFZ54wgnAuq+QlcDC\nhd5n2BrbSy+lX//z5sEDD6Thn2Y9UbhKoF+/tM/w4sV5R2LW+5Ytg4sugj32gP33T+P+nQBsVRSu\nEoCP+gU87M0ayaOPpg1fhg+Hxx+HzTfPOyKrB4WrBMCdw9ZY3nkHvvMd+OpX4bzz4N57nQCsepwE\nzGpURBrrv912acXP6dPhmGPcF2bVVejmILN69eqr8N3vwssvwy23wN575x2R1atCVgKDB3uYqNWn\npUvhwgvTkM/PfS4t+eAEYFkqZCXwqU+loXFm9eT++9N6P1tt5c1erPcUMgmsvz689VbeUZhVx1/+\nAt/7XvrVf9llcOiheUdkjaSQzUFOAlYPliyBMWNgl13SGv/PPOMEYL2vkJXAkCFpzLRZUd13X2r6\n2XbbtN7/ZpvlHZE1qkImAVcCVlQvvQRnnQVPP52afg45JO+IrNFl2hwk6VpJcyRNa/PcepImSpop\n6X5Jg7t73iFD4M03qxurWZYWLoTvfx923z2N/HnmGScAqw1Z9wlcBxzU7rlzgIkRsRUwqfy4W1wJ\nWFGsWAHjxsHWW6cfLs88A//6rzBgQN6RmSWKiGw/QBoB3B0RO5YfPwfsGxFzJG0IlCJimw7eFyuL\n7f3301yBDz7w7EmrXZMnwxlnpMmNl16aKgCzrEkiIiq+MubRJzA0IuaU788Bhnb3BAMHppVEFy2C\ntdaqbnBmq+rFF+Hss9OQz0sugSOP9I8Vq125dgxHREhaaSnS3Nz84f2mpiaampo+fLz++qm8dhKw\nWrFgQRryed11qfP3ppvc7GPZK5VKlEqlHr8/r+agpoiYLWkYMLm7zUGQ1lO/7LI0td4sT++/D//1\nX3DxxXDEEXD++bDhhnlHZY2qu81BeUwWuwsYXb4/GpjQk5NstBG8/nrVYjLrthUr4Be/SJ2+jzwC\nDz0EV13lBGDFkmlzkKSbgX2B9SW9BvwbcCHwa0nfBF4BjurJuTfaKG2wbdbbItKa/uecA2uvDb/6\nFYwalXdUZj2TaRKIiGNX8tIBq3puJwHLw+OPw7/8C8ydm1b7/PKX3elrxVbItYMgJYFZs/KOwhrF\njBnw9a+n3b2OPz7N+D3sMCcAK75CJwFXApa1559PF/399kuzfWfOhH/4B+hbyAVXzD7JScCsA6+8\nAt/8Zhp9ts028MILqRlo0KC8IzOrrsIngYxHuFqD+etf06buu+4Kw4enSuAHP/B8FKtfhU0Ca60F\n/ft7hzGrjtmz4cwz07r+a64Jf/4z/PjHsO66eUdmlq3CJgFI2++9+GLeUViRvfYanHYabLcdtLSk\nDuCLL04z0s0agZOANaSXXoJvfQs++9lUUc6YkWage6KXNRonAWsozz0HJ56Ylh0ZOjSN9rnkEl/8\nrXE5CVhDePppOPpo2Gcf2GqrNNrnxz92s49Z4ZPACy/kHYXVqoi0pv/BB8Pf/m0a5//SS2m0z+Bu\n72dnVp8KPeVlq61SOW/W1vLlcPvtqZln0aK0rPMdd3hZZ7OOFDoJDBuW/oefOxc22CDvaCxv770H\n116bdvEaPhx++MO0ts9qha53zbJV6P89JNhhh7RvqzWuuXPh3/4NRoxIzT833ggPPwxf+YoTgFlX\nCv+/yA47wLRpeUdheXjiCTjppLSe/9y58Pvfp2YfbzRkVrnCJ4Edd3QSaCTLlsGvfw1775128dp2\n2zQ44MorUx+RmXVPofsEAEaOhHHj8o7Csvbmm2nXriuuSKPCzjwzNfd4NU+zVZP5HsMr/WDpFWAh\nsAJYFhF7tHu90z2GWy1ZAuutB3PmpDVfrH5EwJQpcPnlMGECHHkknHpqmuVrZh3r7h7Def6OCtKG\n8/NX5ST9+6eLwh//mNZ8t+J791246Sb4n/+BBQvS8g7PP++JXWZZyLtPoCr7Mo0alTb6tmJ74ol0\nwd90U7jvPrjggjQj/NxznQDMspJ3JfCApBXA/0REj1v2P/e5ND7ciufdd+Hmm9Ov/nnz4OST02Ju\nw4blHZlZY8gzCewVEW9IGgJMlPRcRDzU9oDm5uYP7zc1NdHU1NThiT7/+bQL1NKl0K9fhhFbVbS0\nwEMPwfjxcOed0NQE//EfcOCB0KdP3tGZFUupVKJUKvX4/bl1DH8sCOlHwKKI+Gmb5yrqGG61557p\nQrL//llEaNXw0ktw/fXptsYaMHo0HHecf/WbVVN3O4Zz6ROQNEjSWuX7awBfBFZptP+hh8JvflON\n6Kya3n0Xrrsu/dr/m7+B+fPhttvSqp5nneUEYJa3XCoBSZsBd5Yf9gVujIgL2h3TrUpg6lQ46qi0\noJyq0t1sPbV0aerY/dWv4J57YN9908zeQw5xc51Z1rpbCdREc1BHupsEItLyATfckH5xWu9asQJK\npdTJe+edabvGY4+Fr33Ni/uZ9aaGTQIAY8bAX/6SlhCw7LW0wKOPpgv/rbfCxhunC/9RR8Emm+Qd\nnVljaugkMGsW7LRT2jx8jTUyCqzBLV+eRvbceWeaxbvmmunCf8wx8JnP5B2dmTV0EoDU/LD33nDG\nGRkE1aA++AAmTkwX/rvvTpO5jjgi3bbbzn0wZrWk4ZPAn/4Ehx2WZpr2759BYA1i/nz43e/Shf/+\n+2HnndNF//DD4dOfzjs6M1uZhk8CkFaXHDUKvv/9KgdVxyLgqafgt79Nt6efTqN6Dj88JdUhQ/KO\n0Mwq4SRAqgL22CMNG9100yoHVkfeeQceeCBd9O+9N7XvH3wwfOlLKQF4T16z4nESKLvggnRxmzzZ\na863WrIkjeZ58EGYNCn98t9rr48u/O7YNSs+J4GylpZ0YdtpJ7jkkioGViArVqRqaNKkdOF/5BHY\nZpu0tMb++6cEMGhQ3lGaWTU5CbTx1luwzz5pjZpG6B94//20r8LDD6f9dv/wB9hww48u+vvuC+uu\nm3eUZpYlJ4F2Xn89rVtz2GFw0UX1tUrlnDnp1/3vf58u/NOmwQ47pF/4rbcNN8w7SjPrTU4CHZg3\nL81iXbECrr4attyyKqftVW++mTZdmTLlo7+LFqXVU/feO13w99jDzTtmjc5JYCVWrIDLLktLS5x0\nEpx9NgwdWrXTV82yZfDCC/DMM+k2bVq66L/zDuyyC+y2G+y6a/q7+eaeqGVmH+ck0IVZs1Kz0C9/\nmUbFHH982pu4t4dDvv12Gsr6wgvpNn16uug//3xag2eHHWD77dPfXXeFLbaA1fLeDNTMap6TQIXe\neistdXzTTenX9qhRafXRHXdMt0037XnTygcfwBtvpP6I11//6P6rr6YL/4svpuWWt9gi3bbc8qML\n/jbbuEnHzHrOSaAHFixIyyD/6U8pITzzTKoYBg5Mm56svXa6MA8alJaiWLEi3ZYvT803Cxem5pq3\n305/I9L7hg//6O/w4WllzdYL/5Ahbsoxs+orTBKQdBDwM6APcHVEXNTu9V5LAh2JSOvnvPFG2h1r\n8WJ477004apv3zTKqG9fWH31lCTWWSfdBg9OTUu+wJtZHgqRBCT1Af4MHAD8FfgjcGxEPNvmmFyT\nQNZKpRJNTU15h5GJev5u4O9XdPX+/QqxxzCwB/BCRLwSEcuAXwFfySmWXJRKpbxDyEw9fzfw9yu6\nev9+3ZVXEtgIeK3N41nl58zMrBfllQTqt53HzKxA8uoT2BNojoiDyo/PBVradg5LcqIwM+uBInQM\n9yV1DO8PvA48TruOYTMzy14uK+1HxHJJ/wTcRxoieo0TgJlZ76vZyWJmZpa9mluNRtJBkp6T9Lyk\nutoFQNImkiZLmi7pGUmn5R1TFiT1kTRV0t15x1JtkgZLuk3Ss5JmlPu36oKkc8v/bU6TdJOk/nnH\ntCokXStpjqRpbZ5bT9JESTMl3S9pcJ4xroqVfL9Lyv9tPiXpDknrdHWemkoC5Ulk/wUcBGwHHCtp\n23yjqqplwJkRsT2wJ/DdOvt+rU4HZlCfo8AuA34bEdsCOwF10YwpaQRwMrBLROxIaqY9Js+YquA6\n0rWkrXOAiRGxFTCp/LioOvp+9wPbR8RngZnAuV2dpKaSAHU+iSwiZkfEk+X7i0gXkOH5RlVdkjYG\nDgauBupq8Yzyr6rPR8S1kPq2IuKdnMOqloWkHymDygM3BpFm8xdWRDwELGj39GHA+PL98cDhvRpU\nFXX0/SJiYkS0lB8+Bmzc1XlqLQk0zCSy8i+vnUn/ourJpcDZQEtXBxbQZsCbkq6T9CdJ4yTVxZqv\nETEf+CnwF9KIvbcj4oF8o8rE0IiYU74/B6jBXUWq5hvAb7s6qNaSQD02H3yCpDWB24DTyxVBXZB0\nKDA3IqZSZ1VAWV9gF+DyiNgFeI9iNyd8SNIWwBnACFJ1uqak43INKmPlxcnq8poj6V+BpRFxU1fH\n1loS+CuwSZvHm5CqgbohaXXgduCXETEh73iqbBRwmKSXgZuBL0i6PueYqmkWMCsi/lh+fBspKdSD\n3YBHImJeRCwH7iD9+6w3cyRtCCBpGDA353iqTtJJpCbZipJ4rSWBKcBnJI2Q1A84Grgr55iqRpKA\na4AZEfGzvOOptog4LyI2iYjNSJ2KD0bEiXnHVS0RMRt4TdJW5acOAKbnGFI1PQfsKWlg+b/TA0id\n+/XmLmB0+f5ooK5+iJWX6D8b+EpEfFDJe2oqCZR/gbROIpsB3FJnk8j2Ao4H9isPoZxa/pdWr+qx\n1D4VuFHSU6TRQWNyjqcqIuIp4HrSD7Gny09flV9Eq07SzcAjwNaSXpP098CFwIGSZgJfKD8upA6+\n3zeAnwOJHkCuAAACpklEQVRrAhPL15fLuzyPJ4uZmTWumqoEzMysdzkJmJk1MCcBM7MG5iRgZtbA\nnATMzBqYk4CZWQNzEjAza2BOAmZmDcxJwOqKpAclfbHdc2d0NnNSUuaL+Ek6rbwJzQ1Zf5ZZdzgJ\nWL25mU9uhnI00Nlqir0xbf4fgQMi4oRe+CyzijkJWL25HTikvDFK674NwyPiYUkTJE0pb+15cvs3\nlhcubLtV31mSftTm8fGSHiuvyXKlpE/8/yPpe+XtGadJOr383JXA5sDvJJ1R7S9stir65h2AWTVF\nxHxJj5OW0r2LVBXcUn757yNigaSBwOOSbouI9jtPfex0rXfK24AeBYyKiBXl5qXjgBvaHLMrcBJp\nh7zVgMcklSLiFEl/CzSVN2/5GElfAtYnLZ1+J7A4Il7t4T8Cs25xJWD1qG2T0NHlxwCnS3oS+APp\ngvuZbpxzf2BXYIqkqaQVKDdrd8zewB0R8X5EvEdak3+fzk4qaWtgdETcAFwJnEf97FFgBeBKwOrR\nXcClknYGBkXEVElNpAv5nhHxgaTJwIB271vOx38YDWz3+viIOK+Tzw0+vqOa6Lq/YTRwI3xYxewO\njOviPWZV40rA6k55y87JwHV81CG8NrCgnAC2Afbs4K1zgA0krSepP3AoH13EJwFfkzQEoHzMpu3e\n/xBweHljljVIm5g/1EW4/Uj7+lLer/i9iPi/bnxds1XiSsDq1c2k5pijyo9/B5wiaQbwZ1KTUKsA\niIhlks4HHidtdfrhzloR8aykHwD3lzuElwHfoXwBLx8zVdIvyu8HGFferOXDz+jAONKWnJuUj3lE\n0tci4raefW2z7vGmMmZmDczNQWZmDcxJwMysgTkJmJk1MCcBM7MG5iRgZtbAnATMzBqYk4CZWQNz\nEjAza2D/H6v9nYW8j113AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x104f9df50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "psi = gauss(x,np.sqrt(omega))\n",
    "# Kinetic energy\n",
    "Hpsi = -0.5*d2gauss(x,np.sqrt(omega))\n",
    "# Add potential\n",
    "add_pot_on_basis(Hpsi,VQHO,psi)\n",
    "# Check the exact answer - we won't be able to do this normally !\n",
    "print \"Energy at optimum alpha is: \",integrate_functions(psi,Hpsi,num_x_points,dx)\n",
    "\n",
    "alpha_vals = np.linspace(0.1,10.1,1001)\n",
    "energy = np.zeros(1001)\n",
    "i=0\n",
    "e_min = 1e30\n",
    "alph_min = 0.0\n",
    "for alpha in alpha_vals:\n",
    "    psi = gauss(x,alpha)\n",
    "    norm = integrate_functions(psi,psi,num_x_points,dx)\n",
    "    #if np.abs(norm-1.0)>1e-6:\n",
    "    #    print \"Norm error: \",alpha,norm\n",
    "    Hpsi = -0.5*d2gauss(x,alpha)\n",
    "    add_pot_on_basis(Hpsi,VQHO,psi)\n",
    "    energy[i] = integrate_functions(psi,Hpsi,num_x_points,dx)\n",
    "    if energy[i]<e_min:\n",
    "        e_min = energy[i]\n",
    "        alpha_min = alpha\n",
    "    i=i+1\n",
    "print \"Min E and alph: \",e_min, alpha_min\n",
    "pl.plot(alpha_vals,energy)\n",
    "pl.xlabel(r\"Value of $\\alpha$\")\n",
    "pl.ylabel(\"Energy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the scanning of $\\alpha$ gives us a reasonably accurate answer: the exact answer is $\\sqrt{2}$, and we could find that with a more sophisticated search routine.  The method works, certainly for this simple problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Mixing states\n",
    "\n",
    "The classic examples for the variational approach generally are similar to the example above: we specify a functional form for the trial wavefunction and vary the parameter (or parameters) to find the minimum energy.  But we can also treat the expansion coefficients of a set of basis functions as the variational parameters.  So we might write:\n",
    "\n",
    "$$\n",
    "\\vert \\psi\\rangle = c_1 \\vert \\phi_1 \\rangle + c_2 \\vert \\phi_2\\rangle\n",
    "$$\n",
    "\n",
    "where we are mixing the first two basis states for a system together.  This is what is often done in practical quantum mechanical calculations (though we tend to use several tens, hundreds or thousands of basis coefficients, which requires a more sophisticated approach to minimisation).\n",
    "\n",
    "Notice that we can actually simplify the problem above: we will insist that $\\langle \\psi\\vert \\psi\\rangle = 1$, so we can write $c_1^2 + c_2^2 = 1$.  If we set $c_1 = 1$ and treat $c_2$ as the only parameter, then we can write:\n",
    "\n",
    "$$\n",
    "E = \\frac{\\langle \\psi\\vert \\hat{H} \\vert \\psi\\rangle}{\\langle \\psi\\vert \\psi \\rangle}\n",
    "$$\n",
    "\n",
    "and then solve $\\partial E/\\partial c_2 = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Practical example\n",
    "\n",
    "We will apply an electric field to a perfect square well (so that the potential is just proportional to $x$) and see how well mixing the first two states works as an approximation to the ground state of the new system.  We looked at this in the notebook \"An odd perturbation\"; you should be well used to the square well by now, so I won't explain what we're doing below in great detail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFppJREFUeJzt3XGMnPWd3/H3tyZRFTWcDyGBwQbSFif0FDh8OqCXXDNp\ncsJ2KJyuh1HucqF3RDJQrlFLMFwTHRsR2jMBNT1FdSjJJe4/9UW+6M7nsAGLMroqIjQoa0MOG+NL\nkXaJ7UvSpqZV/vCKb/+YWfZhmLVn55l55pmZ90uymNl9dp4fj3Y/8+zvu9/fLzITSdL0+DujHoAk\nqVoGvyRNGYNfkqaMwS9JU8bgl6QpY/BL0pQpHfwRsTkijkTEyxFx7wrHNCJiLiK+HxHNsueUJPUv\nyvwdf0SsAV4CPgy8CnwX+GhmHi4csxb4NnB9Zi5ExPmZ+eNyw5Yk9avsHf81wLHMfCUzTwN7gJs6\njvkt4M8ycwHA0Jek0Sob/BcD84XnC+2PFV0OnBcRT0fEcxHxOyXPKUkq4ZySX9/LPNHbgE3Ah4B3\nAM9ExHcy8+WS55Yk9aFs8L8KbCg830Drrr9oHvhxZv4M+FlE/BVwFfCm4I8IFw2SpD5kZqzm+LJT\nPc8Bl0fEZRHxduAWYF/HMX8BvD8i1kTEO4BrgRe7vVhm+i+T+++/f+RjqMs/r4XXwmtx5n/9KHXH\nn5mLEXEX8ASwBvhKZh6OiO3tzz+amUci4lvA88DrwGOZ2TX4JUnDV3aqh8ycBWY7PvZox/OHgYfL\nnkuSVJ6duzXUaDRGPYTa8Fos81os81qUU6qBa5AiIusyFkkaFxFBVlzclSSNGYNfkqaMwS9JU8bg\nl6QpY/BL0pQx+CVpyhj8kjRlDH5JmjIGvyRNGYNfkqaMwS9JU6b06pySpOplwlNP9fe1Br8kjZnv\nfQ8+9Sk4fry/r3eqR5LGxPw8fPzj8JGPwC23wAsv9Pc6Br8k1dypU/DpT8Mv/iJccgkcPQrbt8M5\nfc7ZGPySVFOLi/ClL8G73w2vvgqHDsHnPgfvfGe513WOX5JqJhO++U3YsQPWrYPHH4errx7c6xv8\nklQjc3PLhdvPfx62boVY1f5aZ+dUjyTVwPw83HprK+h/8zfh+edbRdxBhz4Y/JI0UsXC7YYN8NJL\ncMcd/Rdue2HwS9IILC7Crl2wcSMsLMDBg63C7bnnDv/czvFLUoWWCrf33LNcuN20qdoxGPySVJGl\njtsTJ+Dhh4dTuO2FUz2SNGRLHbdbt8K2bcMt3PbC4JekIenWcXv77cMt3PbC4JekAess3C513FZR\nuO1F6eCPiM0RcSQiXo6Ie89w3C9HxGJE/EbZc0pSHWXC/v3w3vfC3r0wOwu7d8P69aMe2ZuV+oUj\nItYAXwQ+DLwKfDci9mXm4S7H7QS+BYxoVkuShqfYcTvKwm0vyt7xXwMcy8xXMvM0sAe4qctxvw/s\nBX5U8nySVCvFjts6FG57UTb4LwbmC88X2h97Q0RcTOvNYFf7Q1nynJI0csXC7fr1rY7bMkslV6ls\n8PcS4l8A7svMpDXNU+P3QUk6s25LJT/4YH0Kt70o+970KrCh8HwDrbv+ol8C9kTr957zgS0RcToz\n93W+2MzMzBuPG40GjUaj5PAkaTCKHbcXXTT4pZJ71Ww2aTabpV4jWjfifX5xxDnAS8CHgB8C/wP4\naGdxt3D8V4G/zMxvdPlclhmLJA1LseN2WEsl9ysiyMxVjabUVE9mLgJ3AU8ALwJ/mpmHI2J7RGwv\n89qSNGqde9yOQ+G2F6Xu+AfJO35JdXHqFOzc2ZrLv+OO1k5YdZ3Dr/yOX5ImSbHjtrjHbV1Dv19j\n8IdHkjRcnYXb2dnRFG6rYvBLmmrj1HE7KE71SJpKxY7bm2+enMJtLwx+SVPl1Cn4zGdaUzl1Wiq5\nSga/pKlQ7Lidn2/tcfvAA/DOd456ZNWbovc4SdOoLh23dWLwS5pYnR230zKHfzZO9UiaON06bm+4\nwdBfYvBLmhjFpZI3bBivpZKrZPBLGntLHbfjvFRylXwflDS2lva43bHDwu1qGPySxlKxcDstHbeD\n4lSPpLEyqUslV8nglzQWioXbSy9tddxauO2PwS+p1rrtcTutHbeD4nulpFqy43Z4DH5JtVNcKvmR\nR2DLFufwB8mpHkm10W2pZP9aZ/AMfkkj163jdtqWSq6SwS9pZIqF24WFyd3jtm58P5VUuaXC7Y4d\nsG6dhduqGfySKtW5VLJz+NVzqkdSJZYKt3bcjp7BL2moXCq5fgx+SUPRrePWwm09+J4raaDsuK0/\ng1/SwLhU8nhwqkdSaS6VPF5KB39EbI6IIxHxckTc2+Xzvx0RhyLi+Yj4dkRcWfackurBwu14KhX8\nEbEG+CKwGfhHwEcj4oqOw34A/JPMvBJ4APjPZc4pafQ6O24PHnSP23FS9n35GuBYZr4CEBF7gJuA\nw0sHZOYzheOfBdaXPKekESkWbtetaz3etGnUo9JqlQ3+i4H5wvMF4NozHH8b8HjJc0oaAQu3k6Ns\n8GevB0bEB4HfA9630jEzMzNvPG40GjQajRJDkzQI8/OtefwDB+D+++ETn3AOf5SazSbNZrPUa0Rm\nz9n91i+OuA6YyczN7ed/ALyemTs7jrsS+AawOTOPrfBaWWYskgbr1CnYubM1l3/nna0F1dzusH4i\ngsxc1e9eZf+q5zng8oi4LCLeDtwC7OsY1CW0Qv9jK4W+pPpYXIRdu2DjRve4nVSlfmHLzMWIuAt4\nAlgDfCUzD0fE9vbnHwX+EPh5YFe0JgRPZ+Y15YYtadA6O25nZ+24nVSlpnoGyakeaXTm5uDuu10q\neRyNYqpH0hgr7nFrx+30MPilKdTZcXv0qB2308Tgl6bISnvcWridLr6/S1PAjlsVGfzShCsWbu24\nFTjVI00sC7daicEvTZhi4Xb9epdK1lsZ/NKE6LbHrUslqxvvAaQx5x63Wi2DXxpjLpWsfjjVI42h\n4h6327ZZuNXqGPzSGCkWbi+5pNVxe/vtFm61Oga/NAY6C7cHD9pxq/55nyDV2FLhdseOVsethVsN\ngsEv1VSxcOtSyRokp3qkmil23Fq41TAY/FJNdFsq2cKthsHgl0asuMdtcalkO241LN5LSCPiHrca\nFYNfGgE7bjVKTvVIFSp23LpUskbF4Jcq0K3j1qWSNSoGvzRE3ZZKtuNWo+b9hjQEdtyqzgx+acDm\n5lqF2+PHW4XbLVucw1e9ONUjDUix4/bmm1uFW/9aR3Vk8EslnToFn/mMHbcaHwa/1Kdi4XZ+3qWS\nNT5KB39EbI6IIxHxckTcu8Ixf9z+/KGIsMSlsZYJ+/fDlVfC17/eKuLu3t2625fGQalfRiNiDfBF\n4MPAq8B3I2JfZh4uHLMV+IeZeXlEXAvsAq4rc15pVObm4O67XSpZ463sHf81wLHMfCUzTwN7gJs6\njrkR2A2Qmc8CayPigpLnlSq11HG7ZYtLJWv8lQ3+i4H5wvOF9sfOdsz6kueVKuEet5pEZb99s8fj\nOu+Lun7dzMzMG48bjQaNRqOvQUllLS7CY4/BZz8L11/f6rhd7+2KaqDZbNJsNku9RmT2mt1dvjji\nOmAmMze3n/8B8Hpm7iwc8yWgmZl72s+PAB/IzJMdr5VlxiINwlLhdseO1lLJDz9sx63qLSLIzFVN\nOpa9438OuDwiLgN+CNwCfLTjmH3AXcCe9hvFTztDX6oDl0rWtCg1x5+Zi7RC/QngReBPM/NwRGyP\niO3tYx4HfhARx4BHgTtLjlkaqKWOW5dK1rQoNdUzSE71qGqnTsHOna0mrDvugHvvtflK46efqR47\ndzV1XCpZ084/StPU6Nzj1qWSNa0Mfk0FC7fSMqd6NNHc41Z6K4NfE6mz4/all9zjVlpi8GuiLC7C\nrl2wceObC7fnnjvqkUn14f2PJkJn4XZ21sKttBKDX2PPwq20Ok71aGwVC7culSz1zuDX2HntteXC\n7aWXulSytFoGv8bGUsftxo2wsNAq3D7wgB230mp5j6Tay2x12d5zD6xbZ8etVJbBr1qbm2sVbo8f\nh4cecg5fGgSnelRLS0slb926XLi94QZDXxoEg1+1Uuy43bDBjltpGAx+1UKx43apcGvHrTQc3kdp\npIp73C4VbjdtGvWopMlm8Gtk7LiVRsOpHlXOpZKl0TL4VRmXSpbqweDX0K20x62FW2k0vNfS0LjH\nrVRPBr+GYqnj9sQJeOQR2LLFOXypLpzq0UB1dtweOuRf60h1Y/BrIDo7bo8etXAr1ZXBr1KKhdti\nx61LJUv15f2Y+rJUuC123Fq4lcaDwa9VK3bcfv7zzuFL46bUVE9EnBcRByLiaEQ8GRFruxyzISKe\njoi/jojvR8S/KnNOjc5S4daOW2m8lZ3jvw84kJkbgafazzudBv51Zv4CcB3wLyPiipLnVYVcKlma\nLGWD/0Zgd/vxbuDXOw/IzBOZebD9+P8Ch4GLSp5XFbDjVppMZe/ZLsjMk+3HJ4ELznRwRFwGXA08\nW/K8GiI7bqXJdtbgj4gDwIVdPvXp4pPMzIjIM7zO3wP2Ap9s3/m/xczMzBuPG40GjUbjbMPTgBX3\nuHWpZKl+ms0mzWaz1GtE5opZffYvjjgCNDLzRESsA57OzPd0Oe5twH5gNjO/sMJrZZmxqJyFhdY8\n/pNPwswM3Habc/jSOIgIMnNVt2dl5/j3Abe2H98K/HmXQQXwFeDFlUJfo/Paa63Av+oqO26laVE2\n+P8I+LWIOAr80/ZzIuKiiPhm+5j3AR8DPhgRc+1/m0ueVyUtFW43bnxz4daOW2nylZrqGSSneqpR\nLNyuW9daOdPCrTS++pnq8Rf6KeIet5LARdqmQnGP223b7LiVpp3BP8G67XF7++0WbqVpZ/BPoMVF\n2LVruXB78KAdt5KWee83QTo7bmdnLdxKeiuDf0J0LpXsHL6klTjVM+aKhdulpZJvuMHQl7Qyg39M\nFQu3l15qx62k3hn8Y6bbUskPPGDHraTeeX84JlwqWdKgGPxjYG4O7r7bjltJg+FUT40t7XG7dat7\n3EoaHIO/hoqF2/Xr3eNW0mAZ/DXSrXD74IN23EoaLO8ha8DCraQqGfwj5lLJkqrmVM+IuFSypFEx\n+CvWuVTy0aMulSypWgZ/RToLt0tLJdtxK6lq3mcOWWarWLu0x62FW0mjZvAPUedSyRZuJdWBUz1D\nUOy4tXArqW4M/gEqFm43bLBwK6meDP4BKO5xOz/vHreS6s170RKKHbdLhdtNm0Y9Kkk6M4O/T0uF\n2+PH7biVNF6c6lmlzj1uX3jBwq2k8WLw96izcOtSyZLGVd/BHxHnRcSBiDgaEU9GxNozHLsmIuYi\n4i/7Pd+oFAu3LpUsaRKUueO/DziQmRuBp9rPV/JJ4EUgS5yvUpmwfz+8972wdy/MzsLXvtbaGEWS\nxllk9pfFEXEE+EBmnoyIC4FmZr6ny3Hrga8BDwL/JjP/2Qqvl/2OZdDm5t5cuN2yxTl8SfUUEWTm\nqhKqzB3/BZl5sv34JHDBCsf9B+Ae4PUS56pEseP25ptbHbf+tY6kSXPG0mREHAAu7PKpTxefZGZG\nxFtu1yPiBuBvM3MuIhplBjpMp07BQw+15vLvvLPVceuqmZIm1RmDPzN/baXPRcTJiLgwM09ExDrg\nb7sc9ivAjRGxFfi7wLkR8V8y8+PdXnNmZuaNx41Gg0ajcfb/gxIWF+HLX4bPfhauv75VuHUOX1Kd\nNZtNms1mqdcoM8f/EPCTzNwZEfcBazNzxQJvRHwA+FQd5viXOm537IALL4RHHnGpZEnjqZ85/jJ/\nhf5HwNcj4jbgFWBbexAXAY9l5ke6fM3Iq7dzc3D33S6VLGl69X3HP2jDvuOfn281YD35JMzMwCc+\nYfOVpPFX9V/1jAX3uJWkN5vY4C/ucbuw0CrculSyJE3g6pzFpZIvusg9biWp00QFf7Fw61LJktTd\nREz1FDtub7nFPW4l6UzGOvi77XHrUsmSdGZjGfzFwu3SUsmf+5zLLEhSL8bq3tjCrSSVNzbBv7TH\nrYVbSSqn9lM9S3vcbt0K27ZZuJWksmob/HbcStJw1C74O/e4PXjQjltJGqRa3T/v379cuJ2dtXAr\nScNQq9U5r7giXSpZklahn9U5axX8p0+nc/iStApjvyyzoS9Jw1er4JckDZ/BL0lTxuCXpClj8EvS\nlDH4JWnKGPySNGUMfkmaMga/JE0Zg1+SpozBL0lTxuCXpClj8EvSlOk7+CPivIg4EBFHI+LJiFi7\nwnFrI2JvRByOiBcj4rr+hytJKqvMHf99wIHM3Ag81X7ezX8EHs/MK4ArgcMlzjkVms3mqIdQG16L\nZV6LZV6LcsoE/43A7vbj3cCvdx4QET8H/Gpm/glAZi5m5v8pcc6p4Df1Mq/FMq/FMq9FOWWC/4LM\nPNl+fBK4oMsx7wJ+FBFfjYjvRcRjEfGOEueUJJV0xuBvz+G/0OXfjcXjsrWNV7etvM4BNgH/KTM3\nAf+PlaeEJEkV6HvrxYg4AjQy80RErAOezsz3dBxzIfBMZr6r/fz9wH2ZeUOX16vHHpCSNGZWu/Vi\nmc0O9wG3Ajvb//3zLoM5ERHzEbExM48CHwb+utuLrXbgkqT+lLnjPw/4OnAJ8AqwLTN/GhEXAY9l\n5kfax10FfBl4O/A3wO9a4JWk0ek7+CVJ46nyzt2I2BwRRyLi5Yi4d4Vj/rj9+UMRcXXVY6zK2a5F\nRPx2+xo8HxHfjogrRzHOYevle6J93C9HxGJE/EaV46tSjz8fjYiYi4jvR0Sz4iFWpoefj/Mj4lsR\ncbB9Lf7FCIZZiYj4k4g4GREvnOGY3nMzMyv7B6wBjgGXAW8DDgJXdByzlVbDF8C1wHeqHGPNrsU/\nBn6u/XjzJF6LXq5D4bj/BuwH/vmoxz3C74m1tOpk69vPzx/1uEd4LWaAf790HYCfAOeMeuxDuh6/\nClwNvLDC51eVm1Xf8V8DHMvMVzLzNLAHuKnjmDcawzLzWWBtRHTrERh3Z70WmflMLtdDngXWVzzG\nKvTyPQHw+8Be4EdVDq5ivVyL3wL+LDMXADLzxxWPsSq9XIvjwLntx+cCP8nMxQrHWJnM/O/A/z7D\nIavKzaqD/2JgvvB8of2xsx0ziYHXy7Uoug14fKgjGo2zXoeIuJjWD/2u9ocmtTDVy/fE5cB5EfF0\nRDwXEb9T2eiq1cu1eAz4hYj4IXAI+GRFY6ujVeVmmT/n7EevP7Cdf9o5iT/oPf8/RcQHgd8D3je8\n4YxML9fhC7T6PzIigrd+f0yKXq7F22g1RX4IeAfwTER8JzNfHurIqtfLtfi3wMHMbETEPwAORMRV\nmfnakMdWVz3nZtXB/yqwofB8A613pjMds779sUnTy7WgXdB9DNicmWf6VW9c9XIdfgnY08p8zge2\nRMTpzNxXzRAr08u1mAd+nJk/A34WEX8FXAVMWvD3ci1+BXgQIDP/JiL+J/Bu4LlKRlgvq8rNqqd6\nngMuj4jLIuLtwC20GsGK9gEfB2gv4fzTXF4TaJKc9VpExCXAN4CPZeaxEYyxCme9Dpn59zPzXdnq\nAN8L3DGBoQ+9/Xz8BfD+iFjTXvfqWuDFisdZhV6uxRFaTaG057PfDfyg0lHWx6pys9I7/sxcjIi7\ngCdoVe2/kpmHI2J7+/OPZubjEbE1Io7RWtvnd6scY1V6uRbAHwI/D+xq3+2ezsxrRjXmYejxOkyF\nHn8+jkTEt4DngddpNUtOXPD3+H3x74CvRsQhWjexOzLzf41s0EMUEf8V+ABwfkTMA/fTmvbrKzdt\n4JKkKePWi5I0ZQx+SZoyBr8kTRmDX5KmjMEvSVPG4JekKWPwS9KUMfglacr8f9pB/XNxLB43AAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1089f6d90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define the eigenbasis - normalisation needed elsewhere\n",
    "def eigenbasis_sw(n,width,norm,x):\n",
    "    \"\"\"The eigenbasis for a square well, running from 0 to a, sin(n pi x/a)\"\"\"\n",
    "    fac = np.pi*n/width\n",
    "    return norm*np.sin(fac*x)\n",
    "\n",
    "# We will also define the second derivative for kinetic energy (KE)\n",
    "def d2eigenbasis_sw(n,width,norm,x):\n",
    "    \"\"\"The eigenbasis for a square well, running from 0 to a, sin(n pi x/a)\"\"\"\n",
    "    fac = np.pi*n/width\n",
    "    return -fac*fac*norm*np.sin(fac*x)\n",
    "\n",
    "# Define the x-axis\n",
    "width = 1.0\n",
    "num_x_points = 101\n",
    "x = np.linspace(0.0,width,num_x_points)\n",
    "dx = width/(num_x_points - 1)\n",
    "# Now set up the array of basis functions - specify the size of the basis\n",
    "num_basis = 10\n",
    "# These arrays will each hold an array of functions\n",
    "basis_array = np.zeros((num_basis,num_x_points))\n",
    "d2basis_array = np.zeros((num_basis,num_x_points))\n",
    "\n",
    "# Loop over first num_basis basis states, normalise and create an array\n",
    "# NB the basis_array will start from 0\n",
    "for i in range(num_basis):\n",
    "    n = i+1\n",
    "    # Calculate A = <phi_n|phi_n>\n",
    "    integral = integrate_functions(eigenbasis_sw(n,width,1.0,x),eigenbasis_sw(n,width,1.0,x),num_x_points,dx)\n",
    "    # Use 1/sqrt{A} as normalisation constant\n",
    "    normalisation = 1.0/np.sqrt(integral)\n",
    "    basis_array[i,:] = eigenbasis_sw(n,width,normalisation,x)\n",
    "    d2basis_array[i,:] = d2eigenbasis_sw(n,width,normalisation,x)\n",
    "\n",
    "# Define the potential in the square well\n",
    "def square_well_linear_potential(x,V,a):\n",
    "    \"\"\"Potential for a particle in a square well, expecting two arrays: x, V(x), and potential height, a\"\"\"\n",
    "    for i in range(x.size):\n",
    "        V[i] = a*(x[i]-width/2.0)\n",
    "    # Plot to ensure that we know what we're getting\n",
    "    pl.plot(x,V)\n",
    "    \n",
    "# Declare space for this potential (Vdiag) and call the routine\n",
    "Vdiag = np.linspace(0.0,width,num_x_points)\n",
    "square_well_linear_potential(x,Vdiag,1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now build the Hamiltonian - we could actually evaluate the expectation value by using vector-matrix-vector algebra, though we will use integration below for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Full Hamiltonian\n",
      "   4.935   -0.180   -0.000   -0.014    0.000   -0.004    0.000   -0.002    0.000   -0.001\n",
      "  -0.180   19.739   -0.195   -0.000   -0.018    0.000   -0.006   -0.000   -0.002    0.000\n",
      "   0.000   -0.195   44.413   -0.199    0.000   -0.020   -0.000   -0.006    0.000   -0.003\n",
      "  -0.014   -0.000   -0.199   78.957   -0.200   -0.000   -0.021   -0.000   -0.007   -0.000\n",
      "  -0.000   -0.018   -0.000   -0.200  123.370   -0.201    0.000   -0.021   -0.000   -0.007\n",
      "  -0.004    0.000   -0.020   -0.000   -0.201  177.653   -0.201   -0.000   -0.022    0.000\n",
      "  -0.000   -0.006   -0.000   -0.021    0.000   -0.201  241.805   -0.202    0.000   -0.022\n",
      "  -0.002    0.000   -0.006   -0.000   -0.021   -0.000   -0.202  315.827   -0.202    0.000\n",
      "   0.000   -0.002    0.000   -0.007   -0.000   -0.022    0.000   -0.202  399.719   -0.202\n",
      "  -0.001    0.000   -0.003   -0.000   -0.007    0.000   -0.022    0.000   -0.202  493.480\n"
     ]
    }
   ],
   "source": [
    "# Declare space for the matrix elements\n",
    "Hmat2 = np.eye(num_basis)\n",
    "\n",
    "# Loop over basis functions phi_n (the bra in the matrix element)\n",
    "print \"Full Hamiltonian\"\n",
    "for n in range(num_basis):\n",
    "    # Loop over basis functions phi_m (the ket in the matrix element)\n",
    "    for m in range(num_basis):\n",
    "        # Act with H on phi_m and store in H_phi_m\n",
    "        H_phi_m = -0.5*d2basis_array[m] \n",
    "        add_pot_on_basis(H_phi_m,Vdiag,basis_array[m])\n",
    "        # Create matrix element by integrating\n",
    "        Hmat2[m,n] = integrate_functions(basis_array[n],H_phi_m,num_x_points,dx)\n",
    "        # The comma at the end prints without a new line; the %8.3f formats the number\n",
    "        print \"%8.3f\" % Hmat2[m,n],\n",
    "    # This print puts in a new line when we have finished looping over m\n",
    "    print\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we set up the crude parameter scan.  I have chosen a range for $\\alpha$ which is helpful, though this approach sometimes requires a degree of trial and error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Min E and alph:  4.93261131357 0.012\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x108a0b190>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAERCAYAAACD9ivUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYVPXZ//H3LSuC2MCgCBbEEAFDFFFEsawoCkgUSzRG\nDPj4GExsiTWWKBHFFks0kaDYNRobRCIiENkIFhBdigUUBX8RULGgwiPKwv374zvIsm6Zs1POmZnP\n67r22plT5txzGPaebzd3R0REJIqN4g5AREQKj5KHiIhEpuQhIiKRKXmIiEhkSh4iIhKZkoeIiESW\n8+RhZovMbI6ZVZrZjFr2dzKzl8xslZmdF+VcERGJR1keruFAubt/Vsf+T4GzgIGNOFdERGKQr2or\nq2uHuy9z95nA6qjniohIPPKRPByYbGYzzey0PJ4rIiI5ko9qq17uvtTMWgOTzGyeu0/Nw7kiIpIj\nOU8e7r409XuZmY0BegBpJYCGzjUzTcwlItII7p5Rk0BOq63MbFMz2zz1uAVwGDC3rsMbc66768ed\nK664IvYYkvKje6F7oXtR/0825LrksS0wxszWXeshd59oZkMB3H2UmbUBXgG2ANaa2TlAF2Ab4Mma\n5+Y4XhERSUNOk4e7LwT2qGX7qGqPPwR2qOX0FbWdKyIi8dMI8yJRXl4edwiJoXuxnu7FeroX2WXZ\nqv+Kg5l5IccvIhIHM8OT3GAuIiLFSclDREQiU/IQEZHIlDxERCQyJQ8REYlMyUNERCJT8hARkciU\nPEREJDIlDxERiUzJQ0REIlPyEBGRyJQ8REQkMiUPERGJTMlDREQiU/IQEZHIlDxERCQyJQ8REYlM\nyUNERCJT8hARkciUPEREJDIlDxERiUzJQ0REIlPyEBGRyAo+ebzxRtwRiIiUnoJPHtdfH3cEIiKl\nx9w97hgazcy8VSunshJ23DHuaERECoOZ4e6WyWsUfMnjlFPg5pvjjkJEpLQUfMnjgw+crl1hwQJo\n1SruiEREkk8lD6BdOxg4EP7yl7gjEREpHQVf8nB35s2DAw+EhQuhRYu4oxIRSTaVPFI6dYIDDoA7\n74w7EhGR0lAUJQ+AmTPh6KPh3XehadOYAxMRSTCVPKrZay/o3BkefDDuSEREil/RlDwApkyB00+H\nN9+EJk1iDExEJMEKouRhZovMbI6ZVZrZjFr2dzKzl8xslZmdV2NfXzObZ2bvmNlFDV2rvDx0133y\nySy+ARER+Z6clzzMbCHQ3d0/q2N/a2AnYCDwubvfmNreBJgPHAosBl4BTnT3t6qd6zXjHzcO/vAH\nqKwEyyiviogUp4IoeaTUGaS7L3P3mcDqGrt6AAvcfZG7rwYeAY5q6EIDBoTf48c3OlYREWlAPpKH\nA5PNbKaZnRbhvHbAf6s9/yC1rV5mcMklcNVVUMDNOSIiiVaWh2v0cvelqeqpSWY2z92npnFeWn/6\nhw0b9t3j8vJyysvLOfZYuPzy0IDeu3fjghYRKRYVFRVUVFRk9TXz2tvKzK4AVqxr16hvn5n1BIa5\ne9/U84uBte5+XbVzvtfmsc6998L998Nzz2X/fYiIFLLEt3mY2aZmtnnqcQvgMGBuXYfXeD4T6Ghm\n7c2sKXAC8FS61z7ppDBdyQsvNCJwERGpV05LHma2MzAm9bQMeMjdrzGzoQDuPsrM2hB6Um0BrAW+\nArq4+woz6wfcAjQB7nL3a2q8fp0lD4A77gjddidMyPY7ExEpXNkoeRTVIMGavvkGOnaExx+HHj3y\nGJiISIIlvtoqbptsAhddBMOHxx2JiEhxKeqSB8CqVbDLLmHw4J575ikwEZEEU8kjDc2awQUXwJVX\nxh2JiEjxKPqSB8DXX4fSx9NPQ7dueQhMRCTBVPJIU/PmcOGFKn2IiGRLSZQ8YH3pY/x42GOPHAcm\nIpJgKnlEoNKHiEj2lEzJA1T6EBEBlTwiW1f6qDaXooiINEJJlTwglD46doSxY8O65yIipUYlj0Zo\n3hwuvhiuuCLuSEREClfJlTxg/ZxXjz4KPXvmIDARkQRTyaORNtkELrssLBglIiLRlWTyABgyBBYs\ngOefjzsSEZHCU7LJo2nT0Ovq0ku11rmISFQlmzwgrDb46afw7LNxRyIiUlhKOnk0aRJGnF92mUof\nIiJRlHTyADjmGFi7FsaMafhYEREJSrKrbk3jx4c1P+bMCaUREZFipq66WdKvH2y9NTz4YNyRiIgU\nBpU8Ul54ITSgz58fxoGIiBQrlTyyqFcv6NoVRo2KOxIRkeRTyaOaOXPgsMPgnXdg882z9rIiIomi\nkkeW/eQncOihcNNNcUciIpJsKnnUsHAh7L03vPkmbLNNVl9aRCQRslHyUPKoxW9/C2vWwG23Zf2l\nRURip+SRo+SxbBl07gzTp4dla0VEionaPHKkdetQ+rj00rgjERFJJpU86rByJfzoR2Hakh49cnIJ\nEZFYqOSRQy1ahCnbL7hAkyaKiNSk5FGPU04JU7aPGxd3JCIiyaLkUY+yMrj+erjwQli9Ou5oRESS\nQ8mjAf36Qbt2MHp03JGIiCSHGszTUFkZksjbb8MWW+T8ciIiOaUG8zzp1i0kj2uuiTsSEZFkyHnJ\nw8wWAV8Ca4DV7v69jq9mdivQD/g/YIi7V6Zzbr5KHgCLF4e5r159Fdq3z8slRURyolBKHg6Uu3u3\nOhJHf+CH7t4R+BUwMt1z86ldOzj7bPj97+OMQkQkGfJVbVVfhjsSuA/A3acDW5nZtmmem1fnnw/T\npsGLL8YdiYhIvPJV8phsZjPN7LRa9rcD/lvt+Qepbemcm1ctWsCIEWHqkrVr445GRCQ++Ugevdy9\nG6FN4wwzO6CWY+oqXeyfxrl5NWgQmGm9cxEpbWW5voC7L039XmZmY4AewNRqhywGdqj2fPvUNtx9\nSQPnMmzYsO8el5eXU15envX3UN1GG8Gf/wzHHgvHHAObbZbTy4mIZKyiooKKioqsvmZOe1uZ2aZA\nE3f/ysxaABOBP7r7xGrH9AfOdPf+ZtYTuMXde6Z5bt56W9V08smw445w9dWxXF5EpNESv56Hme0M\njEk9LQMecvdrzGwogLuPSh33F6AvsBI4xd1fM7MOwJM1z63x+rElj8WLYffdYcYM6NAhlhBERBol\n8ckj1+JMHhAaz195JUzbLiLJNW0a3Huvphlap1DGeRStc8+FuXNh4sSGjxWReFRVwRlnQJ8+cUdS\nXJQ8MtCsGdx8cxg8+O23cUcjIrUZNQq23hqOPz7uSIqLqq0y5A5HHAG9e4dBhCKSHMuWwW67wZQp\n4bcEeWnzMLObgLvc/Y1MLpQLSUgeEGbb3W8/mDMH2raNOxoRWee002DzzeGmm+KOJFnylTxOA4YA\nGwN3Aw+7+xeZXDRbkpI8AC69FBYuhL//Pe5IRATg5ZfDWKx587SUQk157W1lZp0ISeQXwDTgTnef\nksnFM5Wk5LFyJXTpEnp0HHxw3NGIlLaqKujRI1Ql/+IXcUeTPHnrbWVmTYBOQGdgGTAbONfM/pHJ\nxYtJixZwyy2hV4eWrBWJ18iRsNVWcOKJcUdSvNKptroZ+CnwHDDa3WdU2zff3XfNbYj1xpaYkgeE\nxvP+/UPj+QUXxB2NSGn68EPo2hWefx46d447mmTKV5vHKcCj7r6yln1bufvyTALIRNKSB8CCBdCz\nJ7z2Wpi+RETy66STYIcd4Npr444kufKVPLoTpkav7gvgfXevyuTimUpi8gAYPjysODh2bNyRiJSW\nSZPgV7+C118PVclSu3wlj5eB7sCc1KauwBvAlsCv3f3ZTALIRFKTxzffhCVrb7gBjjwy7mhESsOq\nVaG66pZbwtgrqVu+GsyXAHu4e3d37w7sAbwH9AGuz+TixWqTTUKD3VlnwYoVcUcjUhpGjAiTlSpx\n5Ec6JY833H232raZ2Sx33yOnEdYfWyJLHuucfDJssw3ceGPckYgUt7feggMOgNmzoV27ho8vdfmq\ntnoU+BR4hLDi3/FAa2AQMM3d984kgEwkPXksWwY//jE88wzsuWfc0YgUp7Vr4aCDwtxVZ50VdzSF\nIV/JozlwBtArtekF4HZgFdDC3b/KJIBMJD15QBg0eNttMH06lOV83UaR0jNqFNxzD7zwAjRpEnc0\nhSHnycPMyoBJ7p7IMdOFkDzc4ZBDYMCAMIW7iGTPkiWhneO550JjuaQnXyWPfwPHxjmeoy6FkDwA\n3nkH9t03LBy1885xRyNSPI49NkwLNHx43JEUlmwkj3QqUlYCc81sUuoxgLv72ZlcuJR07Bjm2Bk6\nFJ59FiyjfzIRAXjiCXjjDXjoobgjKU3plDyGpB6uO9AIyeO+HMaVlkIpeUCY72qffcLCUUOGxB2N\nSGH77LPQGeWxx6BXr4aPlw3lbVZdM9sU2NHd52VysWwrpOQBUFkJhx8e1v1o0ybuaEQK1+DBsOWW\ncOutcUdSmPIySNDMjgQqgQmp593M7KlMLlqqunWDU08NM+8WUM4TSZQJE8KkhyNGxB1JaUtnhPkw\nYB/gcwB3rwQ65DCmonbFFfDmm6G4LSLRfPFFmLvqjjtgs83ijqa0pZM8VtfS02ptLoIpBc2ahT7p\nZ58NH38cdzQiheW886BfP+jTJ+5IJJ3k8YaZnQSUmVlHM7sNeDHHcRW1nj3D1CVnnhl3JCKF45ln\nYPJk+NOf4o5EIL3kcRawG/AN8DDwJfDbXAZVCq68MjScP/po3JGIJN/y5aG66q67YPPN445GIMIa\n5klUaL2tapo+PUzZPnu2el+J1Gfw4LA+x+23xx1JccjXCPNdgfOB9qwfVOju3juTC2dDoScPgEsv\nDQvXjB2rwYMitRkzJizrPHu2FnjKlnwljznASOA1YE1qs7v7q5lcOBuKIXl88w306BHmvRo8OO5o\nRJLl44/D3FWPP67BgNmUr+TxamoRqMQphuQB4RtVnz5h7quddoo7GpFkcIejj4ZOnbQeebblayXB\ncWZ2hpltZ2at1v1kclHZ0O67h7mvBg+GNWsaPl6kFNxzDyxcCH/8Y9yRSG3SKXksYv28Vt9x99jn\nhy2WkgeEpHHwwaEB/fzz445GJF4LFoSZqKdMCXNYSXblbW6rpCqm5AGwaFFo/5g0KZRGREpRVRXs\nvz+ceCKcc07c0RSnnFZbmdmF1R7/rMY+zSqTA+3bh/XOf/EL+PrruKMRicdVV4VJD7WkbLLVWfIw\ns0p371bzcW3P41JsJQ8IjYQnnQQtW8Jf/xp3NCL5NW0aHHccvPYatG0bdzTFK18N5pJHZjByJIwf\nD09p7mIpIZ9/Hr44jR6txFEIlDwSaMst4cEHw3QMixfHHY1I7rmHz/vAgTBgQNzRSDrqSx4/MbOv\nzOwroOu6x+uep3sBM1tkZnPMrNLMZtRxzK1m9o6ZzTaz6tVjfc1sXmrfRWm/qyLQqxf85jcwaJC6\n70rxGz0a3n4brrsu7kgkXTnvbWVmC4Hu7v5ZHfv7A2e6e38z2wf4s7v3NLMmwHzgUGAx8Apworu/\nVe3comvzqG7NmjB4sLwcLr887mhEcmPOHDjkEJg6NQwIlNwrpDaP+oI8ErgPwN2nA1uZWRugB7DA\n3Re5+2rgEeConEeaIE2ahOqr22+H//wn7mhEsm/FCjj+eLjpJiWOQpOP5OHAZDObaWan1bK/HfDf\nas8/SG1rW8f2ktK2bRhpO2iQFo+S4uIeqmb32y+sbyOFpazhQzLWy92XmllrYJKZzXP3qTWO0Xyy\n9ejXLySPQYPCgjhNmsQdkUjm7r4bXn0VZtTaEipJl/Pk4e5LU7+XmdkYQnVU9eSxGNih2vPtCaWM\njWts3yG1fQPDhg377nF5eTnl5eVZijxZhg8P9cJXX632Dyl8s2bB738Pzz+vadbzoaKigoqKiqy+\nZk4bzM1sU6CJu39lZi2AicAf3X1itWOqN5j3BG5JNZiXERrMDwGWADMosQbzmpYuhe7d4f774dBD\n445GpHGWL4e99gojyX/+87ijKU3ZaDDPdcljW2CMhVWOyoCH3H2imQ0FcPdR7j7ezPqb2QJgJXBK\nal+VmZ0JPAs0Ae6qnjhK0XbbwUMPhelLpk+HHXeMOyKRaNauhSFDQlWsEkdh08SIBeiGG+Cxx0LX\nxk02iTsakfSNGAHjxkFFhT67cdKsuiWaPNzD/D8/+AGMGhV3NCLpmTABTj01NJC3K7l+k8lSSOM8\nJIvMQvfd558PI3NFku6998JiZ488osRRLFTyKGDz58MBB8A//xkWzhFJohUrwufzV7/SNOtJoWqr\nEk8eAE8/Hf5TvvKKZiKV5Fm7NlSxtmwZSsmmEV2JoGor4YgjwijdY46BVavijkZkQ1ddBR9+GKbY\nUeIoLip5FAH3sGRnWRk88ID+k0oyPPEE/O53oYG8TZu4o5HqVPIQYH0D+vz5cO21cUcjAjNnwumn\nh/Y4JY7ilI+5rSQPmjcP/1H32SfMTnr00XFHJKVq8eKwqNMdd0C32BerllxR8igibdvC2LHQty9s\nvz3svXfcEUmpWbECjjwSzjxTX2CKndo8itA//xka0V98EXbaKe5opFRUVYUSR5s2cOedantLskKY\n20picNRRsHBh6Ik1bRpstVXcEUmxc4ezz4bVq2HkSCWOUqCSR5Fyh3POgblzw7QQmkdIcun668Oq\nl9OmwRZbxB2NNESDBJU86rVmDZxwAmy8cZiNdyP1rZMcuP9++MMf4IUXQlubJJ+Sh5JHg1atgsMO\nC43nN94YdzRSbCZMCHNWTZkCXbrEHY2kS+M8pEHNmoUG9GefDVULItkyY0ZYe/zJJ5U4SpEazEtA\ny5Yheey/P7RqBf/7v3FHJIXu9ddDl9x77oFeveKORuKg5FEi2rWDiRPhoINCMjn22LgjkkK1cGEY\nS3TTTTBgQNzRSFyUPEpIx45hFt6+fWHTTcNSoCJRLF4Mhx4KF18clkOW0qU2jxLTrVsYhT54cFgK\nVCRdH30EhxwCQ4fCGWfEHY3ETcmjBO27Lzz6KBx/fOheKdKQTz4JJY4TT4QLL4w7GkkCJY8SVV4e\nBnUdfXSYxkSkLusSx4ABcPnlcUcjSaHkUcIOOyys/zFwILz0UtzRSBJ98kmoqurbF0aM0LQjsp6S\nR4k7/PAwQvioo1SFJRtatiwkjv794ZprlDhkQ0oeQt++oQpr4MAwUlhkyZLQrXvAAJU4pHZKHgKE\nKqzHHguN6BMmxB2NxOn99+HAA2HQILj6aiUOqZ2Sh3ynvDxMZTJ4cEgkUnrmzQuJ46yz4JJL4o5G\nkkyDBGUD++0XRqL37w/Ll8Npp8UdkeTLK6+EKUeuvTZ8gRCpj5KHfM/uu8N//hOqsj7+OHwDVdVF\ncfv3v8MYjtGjQwIRaYimZJc6LV0aSiA9e8Jtt0GZvmoUpQcegPPPD1WVBx4YdzSSD1rPQ8kj5778\nMkyi2Lw5PPwwtGgRd0SSLe6hJ9Wdd8L48ZpWvZRoPQ/JuS22CJMptmoVvpUuWRJ3RJIN334Lp54K\nTzwRZhhQ4pColDykQU2bhnUbjjsuVGHNmhV3RJKJTz6BPn3g88/h+eehbdu4I5JCpOQhaTEL03Df\neGP4w6OuvIXp9dfDF4B99w2ljs02izsiKVRq85DIKivDhIonnQRXXglNmsQdkaTj8cfh17+Gm28O\nAwCldKnBXMkjNsuWwc9+FhrSH3wQtt467oikLlVVcNllocPDmDGw555xRyRxU4O5xKZ1a5g0CX78\nY+jePQwwk+RZujRMbvjaazBzphKHZE/Ok4eZNTGzSjMbV8u+lmY2xsxmm9l0M9ut2r5FZjYnde6M\nXMcp0W28MdxwQ6gGOeIIuOWW0P1TkmHy5JDYe/eGZ54JCV8kW3JebWVm5wLdgc3d/cga+24AvnT3\n4Wa2K/BXdz80tW8h0N3dP6vntVVtlRDvvRdGKLduHXpm6Q9VfFavhj/8IQz+u+++sJCTSHWJr7Yy\ns+2B/sBooLZAOwNTANx9PtDezKr/2dGkGAWiQweYOjWMF9hjj/BNV/Jv/nzo1Sv0qpo1S4lDcifX\n1VY3AxcAa+vYPxs4BsDMegA7Adun9jkw2cxmmpmm5ysATZvC9deHBvTTTw89e1aujDuq0rB2Ldx6\na0gcQ4bAuHEq/Ulu5Wy2IjMbAHzs7pVmVl7HYdcCfzazSmAuUAmsSe3b392XpEoik8xsnrtPrfkC\nw4YN++5xeXk55eV1XUry5eCDYc4cOPts6No1TH9xyCFxR1W8FiwIsx9/801YTrhjx7gjkqSpqKig\noqIiq6+ZszYPMxsBnAxUAc2ALYAn3P2X9ZyzEOjq7itqbL8CWOHuN9bYrjaPhBs/PpRCDj88lEpa\ntow7ouJRVRU6K1x3HVx6aUjWGnMj6Uh0m4e7X+LuO7j7zsDPgedqJg4z29LMmqYenwb8x91XmNmm\nZrZ5ansL4DBCyUQKTP/+MHdumJG3S5dQpaV8n7kXXgg9qSZOhBkz4He/U+KQ/MrnOA8HMLOhZjY0\nta0LMNfM5gGHA+ektm8LTDWzWcB04F/uPjGPsUoWbbkljBwJY8fCTTeFaq3Zs+OOqjB9+CH8z/+E\n5YIvuSQkjw4d4o5KSpFGmEteVVWFNpBhw2DgQBg+HLbZJu6oku/rr0MV1Y03wimnwOWXhxmPRRoj\n0dVWIrUpKwu9sObNC2uDdO4cEslXX8UdWTJVVcHdd0OnTmGE+PTp8Kc/KXFI/JQ8JBYtW4YqrJkz\n4d134Yc/DKPV1bU3WLMGHnkk9Fa7//4wL9WTT4b7JJIESh4Sq513DiOhJ08O82N16ADXXgvLl8cd\nWTy+/TaMCu/SJSz9e8stMGUK7Ldf3JGJbEjJQxKha1d49FF47jl44w3YZRc47zx4//24I8uP5ctD\nV+YOHULyGDkSpk0LXZxN8yxIAil5SKLstlsoiaxbrXDPPeGoo0KvorV1zVNQwF57LQzw23nnMLBy\n3LiQQHv3VtKQZFNvK0m0lSvh738P38Q//RQGDw4/u+wSd2SNt2xZaMO4777weOjQsJ54mzZxRyal\nQotBKXmUlFmz4N57QzLZaSc44QQ49tjwrT3pPvkkjHN57DF4+WX46U9DEuzdW4P7JP+UPJQ8SlJV\nFVRUwD/+AU89BT/4Qfhj3KdPaFhu3jzuCENvqcrKUN329NNhlP3hh4fVF/v319rhEi8lDyWPkrd2\nbeil9a9/hR5bc+fC3nvDvvtCz56w116w3Xa5bz/44ov14zBeeik0drdtGyaEPOIIOOggaNYstzGI\npEvJQ8lDavjyS3jxxVA19NJL8OqrIXH85Cew665hnMQuu0C7duGPe+vWYUXEhqxdG6qeli6FJUvC\n4lfvvhvWz5g7Fz77DLp1g332CT8HHgjbbpv79yvSGEoeSh7SAPcwH9ScOfD222H68vfeCwlg8eKQ\nEJo1C/NvNW8eEklZWah2qqoK04J88UVouG/VKpRittsudKnt0AF+9KPQzbh9e9hIfRelQCh5KHlI\nhtxDYli+PKyH8e23YRnXsrLwsy6xbLaZGraleCh5KHmIiESmiRFFRCQWSh4iIhKZkoeIiESm5CEi\nIpEpeYiISGRKHiIiEpmSh4iIRKbkISIikSl5iIhIZEoeIiISmZKHiIhEpuQhIiKRKXmIiEhkSh4i\nIhKZkoeIiESm5CEiIpEpeYiISGRKHiIiEpmSh4iIRKbkISIikSl5iIhIZEoeIiISWc6Th5k1MbNK\nMxtXy76WZjbGzGab2XQz263avr5mNs/M3jGzi3Idp4iIpC8fJY9zgDcBr2XfJcBr7r478EvgzxAS\nDvAXoC/QBTjRzDrnIdaCVVFREXcIiaF7sZ7uxXq6F9mV0+RhZtsD/YHRgNVySGdgCoC7zwfam9k2\nQA9ggbsvcvfVwCPAUbmMtdDpP8Z6uhfr6V6sp3uRXbkuedwMXACsrWP/bOAYADPrAewEbA+0A/5b\n7bgPUttERCQBcpY8zGwA8LG7V1J7qQPgWmArM6sEzgQqgTXUXsUlIiIJYe65+TttZiOAk4EqoBmw\nBfCEu/+ynnMWAl2BHwPD3L1vavvFwFp3v67G8UoyIiKN4O51falPS86SxwYXMTsION/df1pj+5bA\n1+7+rZmdBvRy9yFmVgbMBw4BlgAzgBPd/a2cBysiIg0qy+O1HMDMhgK4+yhCT6p7UyWI14FTU/uq\nzOxM4FmgCXCXEoeISHLkpeQhIiLFJfEjzM2slZlNMrO3zWyimW1Vx3F3m9lHZja3MecXggj3otYB\nlmY2zMw+SA3arDSzvvmLPjvSGTxqZrem9s82s25Rzi0kGd6LRWY2J/U5mJG/qHOjoXthZp3M7CUz\nW2Vm50U5t9BkeC/S/1y4e6J/gOuBC1OPLwKureO4A4BuwNzGnF8IP+m8F0I13wKgPbAxMAvonNp3\nBXBu3O8jg/df53urdkx/YHzq8T7Ay+meW0g/mdyL1POFQKu430ce70VrYC/gKuC8KOcW0k8m9yLq\n5yLxJQ/gSOC+1OP7gIG1HeTuU4HPG3t+gUjnvTQ0wDKjHhYxS2fw6Hf3yN2nE7qCt0nz3ELS2Hux\nbbX9hfxZqK7Be+Huy9x9JrA66rkFJpN7sU5an4tCSB7buvtHqccfAdvWd3AOzk+SdN5LQwMsz0pV\nYdxVgFV46QwereuYtmmcW0gyuRcQOrBMNrOZqZ6OhSyTQcXFNiA50/eT9ucin72t6mRmk4A2tey6\ntPoTd/dMxnZken4+ZOFe1Pf+RgJXph4PB24k1cOtQKT7b1cs36jrk+m92N/dl5hZa2CSmc1Lld4L\nUSb/pxP996ARMn0/vdx9aTqfi0QkD3fvU9e+VCN4G3f/0My2Az6O+PKZnp9XWbgXi4Edqj3fgfDt\nA3f/7ngzGw18b6bjhKvzvdVzzPapYzZO49xC0th7sRjA3Zekfi8zszGE6o5CTR7p3ItcnJtEGb0f\nd1+a+t3g56IQqq2eAganHg8Gxub5/CRJ573MBDqaWXszawqckDqPVMJZ52hgbi3nJ1md762apwgz\nNGNmPYHlqaq+dM4tJI2+F2a2qZltntreAjiMwvssVBfl37ZmSawUPxfrbHAvIn8u4u4dkEbvgVbA\nZOBtYCKrlkFkAAADE0lEQVSwVWp7W+Dpasc9TBiN/g2hzu+U+s4vxJ8I96IfYYT+AuDiatvvB+YQ\nJqQcS2hDif19RbwH33tvwFBgaLVj/pLaPxvYs6H7Uqg/jb0XQAdCL5xZhMG5RX8vCFXB/wW+IHSs\n+X/AZqX4uajrXkT9XGiQoIiIRFYI1VYiIpIwSh4iIhKZkoeIiESm5CEiIpEpeYiISGRKHiIiEpmS\nh4iIRKbkISIikSl5SEkys+fM7LAa235rZrfXc86KPMR1tpm9aWYP5PpaIplQ8pBS9TDw8xrbTgD+\nXs85+ZiO4dfAoe5+ch6uJdJoSh5Sqp4AjjCzMgAzaw+0dfdpZjY2tZ7B67WtaZCadG5utefnm9kV\n1Z4PMrPpqaU8/2Zm3/t/Zmbnmtnc1M85qW1/I8wvNMHMfpvtNyySTYmYkl0k39z9s9Qazf0Js47+\nHPhHavcp7v65mTUHZpjZ4+5e2yqV373cugdm1hk4HtjP3dekqsFOAh6odkx3YAhhuuuNgOlmVuHu\np5vZ4UC5u39W8yJm1g/4AWGa7THA/7n7+428BSIZUclDSln1qqsTUs8BzjGzWcBLhD/UHSO85iFA\nd2CmmVUCvYGdaxyzP/Cku3/t7iuBJ4ED63tRM9sVGOzuDwB/Ay4B9owQl0hWqeQhpewp4GYz6wZs\n6u6VZlZOSAA93X2VmU0BmtU4r4oNv3g1r7H/Pne/pJ7rOhuupWA03J4yGHgIvis17Q3c2cA5Ijmj\nkoeULHdfAUwB7mF9Q/kWwOepxNEJ6FnLqR8B25hZKzPbBBjA+j/+/waOSy3jSeqYHWucPxUYaGbN\nU4vuDKThVfyaEtZdwMw2BVa6+/MR3q5IVqnkIaXuYUK10fGp5xOA083sTcKCOi9VO9YB3H21mV0J\nzCAs+/nmdwe4v2VmlwETUw3lq4HfkPrDnzqm0szuTZ0PcKe7z65+jVrcCRxpZjukjnnRzI5z98cb\n97ZFMqPFoEREJDJVW4mISGRKHiIiEpmSh4iIRKbkISIikSl5iIhIZEoeIiISmZKHiIhEpuQhIiKR\n/X8YKI/JCw8bwgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1086e8250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_alpha = 101\n",
    "alpha_vals = np.linspace(-0.1,0.1,n_alpha)\n",
    "energy2 = np.zeros(n_alpha)\n",
    "i=0\n",
    "e_min = 1e30\n",
    "alph_min = 0.0\n",
    "for alpha in alpha_vals:\n",
    "    psi = basis_array[0] + alpha*basis_array[1]\n",
    "    H_psi = -0.5*(d2basis_array[0] + alpha*d2basis_array[1])\n",
    "    add_pot_on_basis(H_psi,Vdiag,psi)\n",
    "    norm = integrate_functions(psi,psi,num_x_points,dx)\n",
    "    #print alpha, norm\n",
    "    #print Hmat2[0,0] + Hmat2[1,0]*alpha + Hmat2[0,1]*alpha + Hmat2[1,1]*alpha*alpha\n",
    "    energy2[i] = integrate_functions(psi,H_psi,num_x_points,dx)/norm\n",
    "    if energy2[i]<e_min:\n",
    "        e_min = energy2[i]\n",
    "        alpha_min = alpha\n",
    "    i=i+1\n",
    "print \"Min E and alph: \",e_min, alpha_min\n",
    "pl.plot(alpha_vals,energy2)\n",
    "pl.xlabel(r\"Value of $\\alpha$\")\n",
    "pl.ylabel(\"Energy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can compare these results to the perturbation theory results in the odd perturbation notebook - they are almost identical.  Now what happens if we make the field stronger ? For perturbation theory, the simple first order corrections failed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Full Hamiltonian\n",
      "   4.935   -3.603    0.000   -0.288    0.000   -0.079    0.000   -0.033    0.000   -0.017\n",
      "  -3.603   19.739   -3.891   -0.000   -0.368   -0.000   -0.112   -0.000   -0.049    0.000\n",
      "   0.000   -3.891   44.413   -3.970   -0.000   -0.400   -0.000   -0.129    0.000   -0.059\n",
      "  -0.288   -0.000   -3.970   78.957   -4.003   -0.000   -0.417   -0.000   -0.138   -0.000\n",
      "   0.000   -0.368   -0.000   -4.003  123.370   -4.019    0.000   -0.426   -0.000   -0.144\n",
      "  -0.079   -0.000   -0.400   -0.000   -4.019  177.653   -4.029   -0.000   -0.432    0.000\n",
      "   0.000   -0.112   -0.000   -0.417    0.000   -4.029  241.805   -4.035    0.000   -0.436\n",
      "  -0.033   -0.000   -0.129   -0.000   -0.426   -0.000   -4.035  315.827   -4.039    0.000\n",
      "   0.000   -0.049    0.000   -0.138   -0.000   -0.432    0.000   -4.039  399.719   -4.042\n",
      "  -0.017   -0.000   -0.059   -0.000   -0.144    0.000   -0.436    0.000   -4.042  493.480\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEMdJREFUeJzt3W2MXNddx/Hfj6R9EUQJlaM43hS5oFAoKqVFuIGCMgg7\nMn1BFRBUPHVt8qICNY2lSjQtUr01wq4tNULIIaqjNuIFJKp4dGz3wUBHrVAbsOLdmqZWY6gluxRD\nW8qTKnmS/HkxM96b8az3zsw99/H7kVba2ZnMOTna/fu3Z//3XEeEAADt9h1VTwAAkB7FHgA6gGIP\nAB1AsQeADqDYA0AHUOwBoAMWLva2P2r7iu1zma+90vZp21+2/Snbty46DgBgfkUk+8cl7Z742kOS\nTkfED0j629FjAEBFXMRFVba3S3oqIl43enxe0j0RccX2Vkn9iPjBhQcCAMwl1Z797RFxZfT5FUm3\nJxoHAJBD8j/QxvBXB85kAIAK3Zzofa/Y3hoR/2b7Dkn/PvkC2/wDAABziAjP+t+kSvbHJS2PPl+W\n9FfTXhQRfERo//79lc+hLh+sBWvBWqx/XL0aWlkJ3XZb6PHHQy++OH9GXjjZ235C0j2Stti+JOn9\nkj4o6WO275d0UdIvLzoOAHTJ2pq0Z490xx3S2bPS0tJi77dwsY+IX9ngqZ2LvjcAdM1gIB06JB09\nKh05Ii0vS5550+Z6qfbsMYNer1f1FGqDtVjHWqzrylpk0/wzz0h33lncexfSZz/XwHZUNTYA1Mlg\nIB08uJ7m9+zZOM3bVszxB1qSPQBUaHJvvsg0n8VBaABQgcFAOnBA2rVLete7pJMn0xV6iWQPAKVL\nuTe/EZI9AJRkMJA+8AFp585y0nwWyR4ASjBO81u3pt2b3wjJHgASmkzzp06VX+glkj0AJFN1ms8i\n2QNAweqS5rNI9gBQoHGa37ZNWl1d/EybopDsAaAAV68O0/yuXdKDD0onTtSn0EskewBY2Nra8MCy\npaXy+uZnRbIHgDll9+b37Rum+ToWeolkDwBzKetMm6KQ7AFgBlVeBbsIkj0A5NS0NJ9FsgeATTQ1\nzWeR7AHgBrJ9801L81kkewCYYpzms33zTS30EskeAK4zuTdfp4uj5kWyB4CRyTR/8mQ7Cr1EsgcA\nSdXcPapMJHsAndaGTps8SPYAOqvJffOzItkD6JyupPkskj2ATllbk/burcfdo8pEsgfQCYOBdODA\nsNPmgQe6keazSPYAWq/tnTZ5kOwBtFYX9+Y3QrIH0Epd6rTJg2QPoFVI89OR7AG0Bml+YyR7AI03\n7rQhzW+MZA+g0Ujz+ZDsATTStBMqKfQbI9kDaJzs3aO62jc/K5I9gMbIXgXbhrtHlYlkD6ARuAp2\nMSR7ALVG33wxSPYAaotOm+KQ7AHUDp02xSPZA6iVyTTflht+V41kD6AWpqV5Cn1xSPYAKkeaT49k\nD6Ay0zptKPRpkOwBVIJOm3IlLfa2L0r6b0kvSBpExI6U4wGov8FAOnhQOnpUOnJkWPDtqmfVfqmT\nfUjqRcQ3E48DoAFI89UpY8+ef7OBjsueacNVsNUoI9n/je0XJH04Ih5LPB6AmuFMm3pIXezfHBFf\ns32bpNO2z0fEZ8dPrqysXHthr9dTr9dLPB0AZWFvvhj9fl/9fn/h93FELD6bPAPZ+yX9b0R8aPQ4\nyhobQLmyaf7YMdJ8kWwrImb+ZzPZnr3tW2x/1+jz75R0r6RzqcYDUD1OqKyvlNs4t0v6Sw9/b7tZ\n0p9ExKcSjgegQnTa1Ftp2zjXDcw2DtAKg4F06NBwb/7wYfbmU5t3G4craAHMjU6b5uBsHAAzo2++\neUj2AGZCmm8mkj2AXLh7VLOR7AFsanV1mOaXljhvvqlI9gA2NE7z994r7dsnnThBoW8qkj2Aqbh7\nVLuQ7AG8BPeCbSeSPYBrxml+2zbSfNuQ7AFcl+bZm28fkj3QcezNdwPJHuioaSdUUujbi2QPdBAn\nVHYPyR7okOzePGfadAvJHugIOm26jWQPtBydNpBI9kCrkeYxRrIHWih73jxpHhLJHmgdzpvHNCR7\noCWm9c1T6DFGsgdagL55bIZkDzQY94JFXiR7oKHYm8csSPZAw5DmMQ+SPdAgpHnMi2QPNACdNlgU\nyR6oOTptUASSPVBTpHkUiWQP1BBpHkUj2QM1wnnzSIVkD9QEJ1QiJZI9UDHOm0cZSPZAhSb35iny\nSIVkD1RgWqcNhR4pkeyBktFpgyqQ7IGS0DePKpHsgRJkO21WV9myQflI9kBCdNqgLkj2QCLZNM8J\nlagayR4oWPa8+XGap9CjaiR7oECcN4+6ItkDBaDTBnVHsgcWRN88moBkD8yJe8GiSUj2wBzYm0fT\nkOyBGZDm0VQkeyAn0jyaLFmyt73b9nnbz9l+T6pxgNTotEEbJEn2tm+SdFTSTklflfSPto9HxJdS\njAekQqcN2iJVst8h6UJEXIyIgaQnJb010VhA4UjzaJtUe/ZLki5lHl+W9KZEYwGFIs2jjVIV+8jz\nopWVlWuf93o99Xq9RNMBNjcYSIcOSUePSocPDwu+XfWs0HX9fl/9fn/h93FErro825vad0taiYjd\no8fvlfRiRBzOvCZSjA3MI5vmjx0jzaO+bCsiZo4hqfbsz0i6y/Z22y+X9DZJxxONBcxt3DfP3jza\nLsk2TkQ8b/udkj4p6SZJH6ETB3Wzuirt3Ts8b569ebRdkm2cXAOzjYOKXL063Jt/5BHpyBFpeZm9\neTTHvNs4XEGLTllbGxb3paVhmucWgegKzsZBJ2TvBbtvH/eCRfeQ7NF63AsWINmjxbJpnnvBoutI\n9milyatg2bJB15Hs0SqTaf7kSQo9IJHs0SKcNw9sjGSPxuOESmBzJHs0GidUAvmQ7NFI0/bmKfTA\nxkj2aJxs3zydNkA+JHs0xrS+eQo9kA/JHo1A3zywGJI9am1apw2FHpgdyR61RacNUBySPWqHvnmg\neCR71AppHkiDZI9a4F6wQFoke1SONA+kR7JHZcZpftcu0jyQGskeleCESqBcJHuUik4boBoke5SG\nvXmgOiR7JEeaB6pHskdSpHmgHkj2SILz5oF6IdmjcJxQCdQPyR6FmZbmKfRAPZDsUQj65oF6I9lj\nIXTaAM1Assfc6LQBmoNkj5ll9+ZJ80AzkOwxk3Ga37aNThugSUj2yGWy0+bECQo90CQke2yKvnmg\n+Uj22NDVq9d32lDogWYi2WMqOm2AdiHZ4yXomwfaiWSPa0jzQHuR7EGaBzqAZN9xpHmgG0j2HUWa\nB7qFZN9BpHmge0j2HcLdo4DuItl3RPZMG86bB7qHZN9y0860odAD3UOybzFOqAQwliTZ216xfdn2\n2dHH7hTjYLrBQDpwYL3ThhMqAaRK9iHp4Yh4ONH7YwN02gCYJuWevRO+NyZMO6GSQg9gLOWe/QO2\n3y7pjKR3R8S3Eo7VaWtr0vLy+t48RR7ApLmLve3TkrZOeep3JT0q6cDo8e9J+pCk+ydfuLKycu3z\nXq+nXq8373Q6aTCQDh6Ujh6VjhwZbt+Y36eAVun3++r3+wu/jyNi8dncaAB7u6SnIuJ1E1+P1GO3\nWXZv/tgx0jzQFbYVETPHulTdOHdkHt4n6VyKcboo2zfP3jyAvFLt2R+2/aMaduV8RdI7Eo3TKfTN\nA5hX8m2cDQdmGye38d78I48M9+aXl9mbB7pq3m0crqCtucm+edI8gHlwNk5NTV4Fe/IkhR7A/Ej2\nNcRVsACKRrKvEe4eBSAVkn1NTHbaUOQBFIlkXzHOmwdQBpJ9hei0AVAWkn0Fpt0LlkIPICWSfcmy\naZ57wQIoC8m+JHTaAKgSyb4E9M0DqBrJPqFpe/MUegBVINknQqcNgDoh2ReMThsAdUSyLxDnzQOo\nK5J9AcYnVGavgqXQA6gTkv2C2JsH0AQk+zlN65un0AOoK5L9HOibB9A0JPsZZDttuAoWQJOQ7HOi\n0wZAk5HsNzHtvHkKPYCmIdnfQDbNc0IlgCYj2U8xrW+eQg+gyUj2EzhvHkAbkexHOG8eQJuR7EXf\nPID263Sy57x5AF3R2WTPmTYAuqRzyZ7z5gF0UaeSPWkeQFd1ItlzQiWArmt9sqfTBgBanOzpmweA\nda1M9qR5AHipViX77Jk2pHkAWNeaZM+ZNgCwscYne/bmAWBzjU727M0DQD6NTPakeQCYTeOSPWke\nAGbXmGRPpw0AzK8RyZ57wQLAYmqd7CdPqOResAAwn9ome06oBIDi1C7Zc948ABRv7mJv+5dsf9H2\nC7bfOPHce20/Z/u87XvzvufamrRjh/T008O9+T17JHveGQIAxhZJ9uck3SfpM9kv2n6tpLdJeq2k\n3ZL+yPYNx+l633y/3696CrXBWqxjLdaxFoubu9hHxPmI+PKUp94q6YmIGETERUkXJO3Y6H2yaf7s\nWWnv3u6leb6R17EW61iLdazF4lLs2W+TdDnz+LKkqbvuXU7zAFCmG3bj2D4taeuUp94XEU/NME5M\n++I4zVPkASAtR0ytw/nfwP60pHdHxDOjxw9JUkR8cPT4E5L2R8TTE//dYgMDQEdFxMyb3UX12WcH\nPi7pT20/rOH2zV2S/mHyP5hnsgCA+SzSenmf7UuS7pZ00vbHJSkinpX0MUnPSvq4pN+ORX99AAAs\nZOFtHABA/SW/gtb27tHFVc/Zfs8Gr/nD0fNrtt+Qek5V2WwtbP/aaA2+YPvvbf9IFfMsQ57vi9Hr\nftz287Z/ocz5lSnnz0jP9lnb/2S7X/IUS5PjZ2SL7U/YXh2txZ4Kppmc7Y/avmL73A1eM1vdjIhk\nH5Ju0rDPfrukl0lalfRDE695i6RTo8/fJOnzKedU1UfOtfgJSd89+nx3l9ci87q/k3RC0i9WPe8K\nvy9ulfRFSXeOHm+pet4VrsWKpEPjdZD0DUk3Vz33BGvx05LeIOncBs/PXDdTJ/sdki5ExMWIGEh6\nUsOLrrJ+XtIfS1IMO3ZutX174nlVYdO1iIjPRcR/jR4+LamtTal5vi8k6QFJfybpP8qcXMnyrMWv\nSvrziLgsSRHx9ZLnWJY8a/E1Sa8Yff4KSd+IiOdLnGMpIuKzkv7zBi+ZuW6mLvZLki5lHk+7wGra\na9pY5PKsRdb9kk4lnVF1Nl0L20sa/qA/OvpSW/+4lOf74i5Jr7T9adtnbP9GabMrV561eEzSD9v+\nV0lrkh4saW51M3PdTH3Ecd4f0Mk2zDb+YOf+f7L9M5J+U9Kb002nUnnW4g8kPRQRYdu6/nukLfKs\nxcskvVHSz0q6RdLnbH8+Ip5LOrPy5VmL90lajYie7e+XdNr26yPifxLPrY5mqpupi/1XJb0q8/hV\neulRCtNec+foa22TZy00+qPsY5J2R8SNfo1rsjxr8WOSnhzWeW2R9HO2BxFxvJwplibPWlyS9PWI\n+Lakb9v+jKTXS2pbsc+zFj8p6fclKSL+2fZXJL1G0plSZlgfM9fN1Ns4ZyTdZXu77ZdreBrm5A/r\ncUlvlyTbd0v6VkRcSTyvKmy6Fra/V9JfSPr1iLhQwRzLsulaRMT3RcSrI+LVGu7b/1YLC72U72fk\nryX9lO2bbN+i4R/kni15nmXIsxbnJe2UpNEe9Wsk/Uups6yHmetm0mQfEc/bfqekT2r4l/aPRMSX\nbL9j9PyHI+KU7bfYviDp/yTtTTmnquRZC0nvl/Q9kh4dJdpBRGx4YmhT5VyLTsj5M3J+dOzIFyS9\nKOmxGF682Co5vy8OSnrc9pqGYfV3IuKblU06EdtPSLpH0pbRxav7NdzOm7tuclEVAHRA7W5LCAAo\nHsUeADqAYg8AHUCxB4AOoNgDQAdQ7AGgAyj2ANABFHsA6ID/B8tmCmxudUz4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10906d2d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Vdiag2 = np.linspace(0.0,width,num_x_points)\n",
    "square_well_linear_potential(x,Vdiag2,20.0) # A much larger potential\n",
    "# Declare space for the matrix elements\n",
    "Hmat3 = np.eye(num_basis)\n",
    "\n",
    "# Loop over basis functions phi_n (the bra in the matrix element)\n",
    "print \"Full Hamiltonian\"\n",
    "for n in range(num_basis):\n",
    "    # Loop over basis functions phi_m (the ket in the matrix element)\n",
    "    for m in range(num_basis):\n",
    "        # Act with H on phi_m and store in H_phi_m\n",
    "        H_phi_m = -0.5*d2basis_array[m] \n",
    "        add_pot_on_basis(H_phi_m,Vdiag2,basis_array[m])\n",
    "        # Create matrix element by integrating\n",
    "        Hmat3[m,n] = integrate_functions(basis_array[n],H_phi_m,num_x_points,dx)\n",
    "        # The comma at the end prints without a new line; the %8.3f formats the number\n",
    "        print \"%8.3f\" % Hmat3[m,n],\n",
    "    # This print puts in a new line when we have finished looping over m\n",
    "    print"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ground state energy:  4.08338931328\n",
      "Ground state wavevector:  [  9.73020341e-01   2.29548881e-01   2.26424844e-02   4.99446267e-03\n",
      "   8.93004688e-04   5.21069181e-04   1.27824770e-04   1.14624515e-04\n",
      "   3.24084063e-05   3.62462861e-05]\n",
      "Min E and alph:  4.10605694331 0.24\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10926b790>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAERCAYAAAB2CKBkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xnc1XP6x/HX1UYlLbJOyJYIbZYYdDDIVmTsa2QGWZKd\nhmb4xRj7UkaUSsreQtNC7pRMhrKVfRlCi1Giwp37+v3xOXK73dW57/uc8znL+/l43I/O8j3ne3Ue\n3/u+zme7PubuiIhI8aoVOwAREYlLiUBEpMgpEYiIFDklAhGRIqdEICJS5JQIRESKXMYSgZkNNrMF\nZvZmhcfPN7O3zewtM/t7ps4vIiKpyWSLYAjQpfwDZrYf0BXYxd13Am7O4PlFRCQFGUsE7j4NWFzh\n4XOAG9y9NHnMokydX0REUpPtMYLtgH3N7N9mVmJmu2b5/CIiUkGdCOdr6u6dzGw34FFg6yzHICIi\n5WQ7EcwDngRw9/+YWZmZbeDu/yt/kJmpAJKISDW4u1X1NdnuGhoN7A9gZq2AehWTwM/at3c++8xx\n109Nf6699troMRTSjz5PfZa5+lNdmZw+OhKYAbQys8/MrAcwGNg6OaV0JHDq6l5//PHQqRO8/HKm\nIhQREchg15C7n7Cap05J5fWXXQatW8Nhh8GAAXDMMWkMTkREVsn2GEGVdO0KkyfDEUfARx+F5GBV\n7v2SRCIRO4SCos8zffRZ5garSb9SppiZl49r3jw4/HDYbbfQOqhbN2JwIiI5yszwPBgsrpYWLWDa\nNPjii5AQvv02dkQiIoUjLxIBQKNGMGYMbLkldO4MX34ZOyIRkcKQN4kAoE4d+Oc/4aijYK+94J13\nYkckIpL/cnqwuDJm8Je/hO6iRAJGjw7TTEVEpHryqkVQXo8e8MADYUbRv/4VOxoRkfyVt4kAwhqD\nsWPh9NNhxIjY0YiI5Ke86xqqaM894fnnoUsX+PprOP/82BGJiOSXvE8EADvuCC+8AAceCEuWQN++\nWngmIpKqvFhQlqr58+Ggg0JCuPlmJQMRKS4FvaAsVZtsAiUlMGMG/OlP8NNPsSMSEcl9BZUIAJo1\ng0mT4P33wyDyypWxIxIRyW0FlwggrEIePx4WLgzlrH/8MXZEIiK5qyATAUCDBmFq6Y8/wtFHww8/\nxI5IRCQ3FWwiAFhnHXjiifBv9+7w/fexIxIRyT0FnQgglKweORIaNgw1ipQMRER+reATAYRk8PDD\n0LgxdOumZCAiUl5RJAIIlUsfeijMKjrqKI0ZiIj8rGgSAYRkMHx4mFWkAWQRkaCgVhanqrQUjjsO\nysrgsce09aWIFAatLK6CunVh1KiQCE48UYvORKS4FWUiAKhXL7QGli6FM84ISUFEpBgVbSKAsL7g\nqafgv/+Fc8+FHOwlExHJuKJOBBBWID/9NMyeDRdfrGQgIsWn6BMBhFlEEybAc8/B9dfHjkZEJLsK\nYmOadGjaFCZOhH32gSZNtNOZiBQPJYJyNtkEJk/+JRmcckrsiEREMi9jXUNmNtjMFpjZm5U8d7GZ\nlZlZs0ydv7patgwtg0svhXHjYkcjIpJ5mRwjGAJ0qfigmW0OHAj8N4PnrpEddwwlrM84A158MXY0\nIiKZlbFE4O7TgMWVPHUrcFmmzpsuu+8eahN17w5z5sSORkQkc7I6a8jMugHz3P2NbJ63ug4+GG69\nFbp0CWsNREQKUdYGi82sAXAVoVto1cPZOn91nXQSLFoEhxwC06eH6qUiIoUkm7OGtgFaAq+bGUAL\n4FUz293dF1Y8uF+/fqtuJxIJEolEVoKsTO/eMG9e2Mtg8mRYd91ooYiIrFJSUkJJSUmN3yej1UfN\nrCUwzt13ruS5j4GO7v51Jc9ltPpodZSVhdZBaSk88gjUrh07IhGRX8u56qNmNhKYAbQys8/MrEeF\nQ3LrL/1a1KoFDz4IX38NffqoFIWIFI6i3I+gJpYsgb33hp49Q5eRiEiuqG6LQCuLq6hJExg/Hvba\nC7bcMmx7KSKSz5QIqmGLLWDMmDCtdLPNYI89YkckIlJ9qj5aTR07wuDBcOSR8PHHsaMREak+JYIa\nOOIIuOoqOPxw+Oab2NGIiFSPBovT4Lzz4P334ZlnoI4620QkkpybPlpMbr8dzODCCzWtVETyjxJB\nGtSpExaZlZTAXXfFjkZEpGrUkZEmjRuHvY/33BNat4aDDoodkYhIatQiSKOttoJHH4WTT4Z3340d\njYhIapQI0mzffaF/f+jaFRZXthuDiEiO0ayhDOndG+bODauQNZNIRLJBs4ZyzM03h38vvzxuHCIi\na6NEkCF16sCoUaEUxfDhsaMREVk9dQ1l2Jw5sN9+oYto111jRyMihUxdQzmqTRu47z7o3h3mz48d\njYjIbykRZMGRR8IZZ8Cxx4YdzkREcom6hrKkrCzsebzVVnDnnbGjEZFCpK6hHFerVhg0njBBg8ci\nklvUIsiyOXMgkYCJE6FDh9jRiEghUYsgT7RpAwMGwNFHw9dfx45GREQtgmguuSS0Dp55JnQbiYjU\nlFoEeebGG2H5cvjb32JHIiLFTi2CiObPD4vMBg2CQw6JHY2I5LvqtgiUCCKbPj2MF8ycCS1bxo5G\nRPKZuoby1N57h8J0xxwDP/wQOxoRKUZqEeQA99Aq2HRTuOee2NGISL5SiyCPmcGQIWFtwciRsaMR\nkWKjFkEOee01OPBAeOEF2GGH2NGISL5Ri6AAtGsHN9wQitMtXx47GhEpFhlNBGY22MwWmNmb5R77\nh5m9bWavm9mTZtY4kzHkmzPPhLZt4YILYkciIsUi0y2CIUCXCo9NAtq4e1vgPeDKDMeQV8xg4ECY\nNk3F6UQkOzKaCNx9GrC4wmOT3b0seXcm0CKTMeSjRo3gscegTx94++3Y0YhIoYs9RnAGMD5yDDlp\nl12gf/8wXrBiRexoRKSQ1Yl1YjO7GvjR3R+u7Pl+/fqtup1IJEgkEtkJLIf07AlTpsBFF8G998aO\nRkRyTUlJCSUlJTV+n4xPHzWzlsA4d9+53GOnA2cBB7j795W8piinj1Zm6dKwb8ENN4TVxyIiq5M3\n00fNrAtwKdCtsiQgv7b++jBqFPTqBR99FDsaESlEGW0RmNlIoDPQHFgAXEuYJVQP+Hlblpfc/dwK\nr1OLoILbb4eHHw5F6urVix2NiOQiVR8tcO7QtSvsuCP8/e+xoxGRXKREUAS++iqsPh48GA46KHY0\nIpJr8maMQKqveXMYOhR69ICFC2NHIyKFQi2CPHTVVTB7tvY7FpFfU4ugiPz1r7B4MdxxR+xIRKQQ\nqEWQpz76CPbYAyZPDuMGIiJqERSZrbeG226DE09UyWoRqRm1CPLcSSeFRWcDB8aORERiU4ugSA0Y\nELa4HDMmdiQikq/UIigAM2ZA9+5hJtGmm8aORkRiUYugiO21F5x9dlhfUFa29uNFRMpTIigQffvC\nN9/A3XfHjkRE8o26hgrIhx9Cp05hD4Odd1778SJSWNQ1JGyzDdx0U5hJ9L0KfItIitQiKDDuYQOb\nli3h5ptjRyMi2aTqo7LKV19B27bw0EOw336xoxGRbFHXkKzSvDk88ACcfjosWRI7GhHJdWoRFLDz\nzgvF6UaMiB2JiGRDxloEZnarmbWpXlgS0003waxZYc9jEZHVSaVr6G3gPjN72czONrPGmQ5K0qNB\nAxg+HC64AD7/PHY0IpKrUu4aMrPWwOnAicB0YJC7P5+RoNQ1lFbXXQfTpsGECdrIRqSQZXSw2Mxq\nA62BHYBFwOtAHzN7pKonlOy78kpYujQUqBMRqWitLQIzuw04ApgC3O/uL5d77l133z7tQalFkHbv\nvRdqEk2fDq1bx45GRDIhY+sIzKwH8Ki7L6vkuSbunvYJikoEmTFwIAweHKqV1q0bOxoRSbdMJoKO\nQMWDvgH+6+4rq3rClIJSIsgIdzjkkNAyuOaa2NGISLplMhH8G+gIvJF8aGdgDtAYOMfdJ1b1pGsN\nSokgYz7/HNq3h/HjYdddY0cjIumUycHiL4B27t7R3TsC7YCPgAOBm6p6Qonrd7+DO+6AU0+FFSti\nRyMiuSCVRLC9u8/5+Y67zwVau/uH/LbLSPLA8ceHMtVXXx07EhHJBakkgjlmNtDMOptZwswGAHPN\nbB2gNMPxSQaYhamkjzwCU6fGjkZEYktljKA+0Av4ffKhF4EBwPdAQ3f/djWvGwwcBix0952TjzUD\nHgG2BD4Bjq1s1pHGCLJj/Hjo1QveeAMaNYodjYjUVEYGi82sDjDZ3atczNjM9gG+A4aVSwQ3AV+5\n+01mdjnQ1N2vqOS1SgRZctZZ4d9Bg+LGISI1l5HB4uT00DIza1LVN3b3acDiCg93BYYmbw8Fjqzq\n+0p63XILPPtsaB2ISHGqk8Ixy4A3zWxy8jaAu/sF1Tjfxu6+IHl7AbBxNd5D0mj99WHIEDj55NBF\n1KxZ7IhEJNtSSQRPJn9+7qsx0jBbyN3dzFb7Pv369Vt1O5FIkEgkanpKWY1EImxv2asXjBwZOxoR\nSVVJSQklJSU1fp+Uqo+aWQNgC3d/p0pvbtYSGFdujOAdIOHu881sU+B5d/9N5RuNEWTfihVhodl1\n14WkICL5J5Mb03QFZgMTkvfbm9nYqocIwFjgtOTt04DR1XwfSbP69WHoUDj/fFiwYO3Hi0jhSGX6\n6Cxgf8K39/bJx95y953W8rqRQGegOWE84BpgDPAosAWaPpqTrr4a3noLRo8O6w1EJH9kstbQTHff\nw8xml0sEb7j7LtWMde1BKRFE88MPsNtucPHFcNppaz9eRHJHJmsNzTGzk4A6Zradmd0FzKhyhJIX\n1lkHhg2DSy+FefNiRyMi2ZBKIjgfaAP8AIwElgK9MxmUxNWuXRgrOPPMULpaRApbynsWZ5O6huIr\nLQ37FvTsCX/+c+xoRCQVmRwj2B64BGjJL+sO3N33r+rJUg5KiSAnzJ0LnTvDyy/DVlvFjkZE1iaT\nieANYCAwC/gp+bC7+6tVjjLVoJQIcsbNN8PTT8OUKVArlY5EEYkmk4ng1eSGNFmjRJA7fvoptAqO\nOQYuvDB2NCKyJplMBP2ARYQyEz/8/Li7f13Vk6UclBJBTnn/fdhzz7DpfatWsaMRkdXJZCL4hEpq\nC7l7xnqNlQhyz113wcMPw/TpULt27GhEpDIZSwQxKBHknrIy+MMf4OCD4fLLY0cjIpVJ+4IyM7us\n3O1jKjzXv6onkvxWqxYMHhwGj996K3Y0IpJOa5oHckK521dVeO6QDMQiOa5lS+jfP5SeKNVu1SIF\nQxMCpUp69oSNNgoJQUQKgxKBVIkZ3H8/3HMPzJoVOxoRSYfVDhab2U/A8uTd+sCKck/Xd/dUdjer\nXlAaLM55w4fDTTfBK6+EQnUiEp9mDUlWuUP37tC6NdxwQ+xoRASUCCSCBQugbduwiU2nTrGjEZFM\n7kcgUqmNN4a77w6ziJYvX/vxIpKb1CKQGjvxxDCT6PbbY0ciUtzUIpBo7r4bHnsMnn8+diQixams\nDPr1q/7rlQikxpo1g/vugzPOgKVLY0cjUlyWLYNjj4VJk6r/HkoEkhaHHQYHHAB9+sSORKR4fPop\n7L03rLdezVrkSgSSNrfeCs89B888EzsSkcL30kthtt5JJ8GQITVbz6PBYkmrqVPD4PEbb8AGG8SO\nRqQwjRgBF10UCkEefvgvj2sdgeSMPn3giy9g1KjYkYgUlrIyuOaasDfI2LGw006/fl6JQHLGihXQ\nsWO4YI8/PnY0IoVh2bKwZmf+fHjqKdhww98eo+mjkjPq14dhw8Iex198ETsakfz3+eew777QsGEY\nh6ssCdSEEoFkxK67wrnnwplnhrpEIlI9s2aFQeFjjoEHH8xMkUclAsmYq66Cr76Cf/4zdiQi+Wn0\n6LA97B13wBVXhDLwmRBljMDMrgROBsqAN4Ee7v5Duec1RlAg3nkH9tknTHXbdtvY0YjkB/ewLewd\nd4RksOuuqb0ub8YIzKwlcBbQwd13BmoDGlIsUK1bw1/+AqeeCitXxo5GJPeVlsKf/hSmiP7736kn\ngZqI0TW0FCgFGphZHaAB8HmEOCRLzjsPGjSAv/89diQiuW3xYujSJcwMmj4dWrTIznmzngjc/Wvg\nFuBT4Atgibs/m+04JHtq1QqDXHfcAa++Gjsakdz00Uew116w886hO2i99bJ37oxtN7k6ZrYN0Bto\nCXwDPGZmJ7n7iPLH9StXSi+RSJBIJLIXpKRdixYhEZx8cpgFUb9+7IhEcsdLL4Ud//r2hV69Un9d\nSUkJJSUlNT5/1geLzew44EB375m8fwrQyd17lTtGg8UF6sQToXlzuPPO2JGI5IZHHw3dpw8+CIce\nWrP3ypvBYuAdoJOZ1TczA/4AzI0Qh0Rwzz2h2TtxYuxIROJyhxtvhEsugcmTa54EaiLW9NHLgNMI\n00dnAT3dvbTc82oRFLApU8IsotdeC60DkWJTWhoWXL76KowbB7/7XXreV7WGJK9ceil88AE8+WTm\nFsmI5KJvvgmrhOvVC4UZ0zkonE9dQyJcfz18/HEooytSLH7eSKZVq+zPDFoTtQgkmjlzIJGAGTNg\nu+1iRyOSWbNmQdeucPHF0Lt3ZlrC6hqSvHTXXTB8OLz4ItStGzsakcx4+mno0SPU3erePXPnUSKQ\nvOQedlhq2xb6948djUj6DRgA110XuoL22COz51IikLy1cCG0axdqq+y3X+xoRNKjrAwuuyy0BsaP\nh623zvw5lQgkr02YEAptvfYaNGsWOxqRmlmxIkyRXrgw7CaWrWtas4Ykr3XpAn/8I5x1ljaykfy2\naBEccEAY85o0KT++2CgRSM644YZQeEsb2Ui+ev/9UDgukYCHHsrMbmKZoK4hySnvvhvmWU+ZEqow\niuSLGTPCjKDrrgst2xjUNSQFYfvtw85Mxx0Hy5fHjkYkNY89BkceGQrHxUoCNaEWgeQc9zDQVr8+\n3Hdf7GhEVs8dbrkllFgfNy7MfotJs4akoHz7LXToEJrZx2sjU8lBK1fChRfCtGnwzDOw+eaxI1Ii\nkAI0axYcfLBKUEjuWbYMTjghTBN9/HFo3Dh2RIHGCKTgdOgAf/0rHHssfP997GhEgvnzoXNn2GCD\nsFAsV5JATSgRSE475xzYdlvo0yd2JCKhUGKnTtCtW6icWyj1sZQIJKeZwf33hx3NHn00djRSzJ5/\nPpRAue46+MtfCmsfDY0RSF74ebxg+vQwxVQkm4YNC1tKjhoF++8fO5rV02CxFLz77gub3s+cCQ0b\nxo5GioE7/O1vYX3AM8/AjjvGjmjNlAik4LnD6aeHqo7DhhVW01xyz48/hkKIc+aENQKbbBI7orXT\nrCEpeGYwcGCoUKqFZpJJixeHQoiLF0NJSX4kgZpQIpC80qABPPEE9O0LL78cOxopRB9/HArH7bIL\nPPlkcXRDKhFI3mnVCgYNCmWrFy6MHY0Ukpkz4fe/h3PPhdtvh9q1Y0eUHRojkLzVt2+YRTR5cuHM\n55Z4nngCzj47rA844ojY0VSPBoul6Pz0U9jvuHVruO222NFIvnKHf/wD7roLxo6F9u1jR1R91U0E\ndTIRjEg21K4NDz8Mu+4KHTvCySfHjkjyTWkp9OoVxpteeglatIgdURxKBJLXmjaFMWPCis/tt4fd\ndosdkeSLxYvhmGNg3XVDBdFGjWJHFI8GiyXv7bRTKEPRvTt8+WXsaCQffPgh7Lln2AVvzJjiTgKg\nRCAFols3+POf4aijVKlU1uyFF8LMoN69w9hSscwMWpMog8Vm1gS4H2gDOHCGu/+73PMaLJYqcw9b\nXK6zjlYeS+WGDIHLLw8byx90UOxo0i+vZg2Z2VBgqrsPNrM6QEN3/6bc80oEUi3Ll0MiAV27huml\nIhDKklx5ZZgiOm4c7LBD7IgyI29mDZlZY2Afdz8NwN1XAt+s+VUiqWnQIPT5duoU9jHQNpfy7bdh\nRtmSJWHB2AYbxI4o98QYI9gKWGRmQ8xslpkNMrMGEeKQArXppmE++Pnnh20upXh98kkYD9hoo7Dw\nUEmgcjESQR2gAzDA3TsAy4ArIsQhBaxt21A6+Oij4YMPYkcjMUybFmYG9ewZihTWqxc7otwVYx3B\nPGCeu/8nef9xKkkE/fr1W3U7kUiQSCSyEZsUkMMOg2uvhUMOCS2DDTeMHZFky333hTGi4cPDhkaF\nqqSkhJKSkhq/T6zB4heAnu7+npn1A+q7++XlntdgsaTN1VfDs8/ClCnFUUmymJWWhmmhU6aEsaJW\nrWJHlF35NmuoLWH6aD3gQ6CHZg1Jpvy8oc3//gejR0MdracvSIsWhZXC660HI0ZA48axI8q+vEoE\na6NEIOlWWhoqSm62WViFXEtLKQvKq6+GleUnnRQ2ly/WRWLaoUxkDerWhccfh7ffDpuQ63tG4Rg2\nLOwmduut0L9/8SaBmlCLQIrK4sXQuTMce6wWnOW7H38MSf1f/wpdfm3axI4ovrxZUCYSU9OmMGkS\n7L03NGkC550XOyKpji++COMBzZqFEtJNm8aOKL+pa0iKziabhFlE//hHmGYo+WXq1LAHxaGHhplB\nSgI1pxaBFKWWLeG558I+BnXqwBlnxI5I1qasDG66KewlPGxYYRaNi0WJQIrWttuGZLD//iEZnHpq\n7Ihkdb7+Gk47LUwB/s9/YPPNY0dUWNQ1JEWtVavQTXTllaFEseSemTOhQwfYbjsoKVESyAS1CKTo\ntW4dVqIeeGAoY92rV+yIBEJX0C23wM03w733hk2HJDOUCEQI+x1PnQoHHBCSwaWXxo6ouC1cGLqC\nli4Ns4K23DJ2RIVNXUMiSVttFbYxfOCBUJ9IS1nimDAB2rWD9u1DV5CSQOZpQZlIBYsWhcqlO+8c\nuiTq1o0dUXH4/vuwjeRTT4VZQSo4XHUqMSGSJhtuCM8/D/Pnw5FHwrJlsSMqfLNnw267wZdfwuuv\nKwlkmxKBSCUaNgxlCzbeOKw1+PLL2BEVppUr4frrw54Bl10GjzyiBWIxKBGIrEbdumG8oFs32GMP\nmDUrdkSF5e23wzaSL7wQqoeecgpYlTs1JB2UCETWwCwMHN92W/jW+vjjsSPKf6Wl8H//B/vsE/aJ\nmDhRawNi0/RRkRQcfTRsvXUYM3jlldCdoQ1uqm7WLDjzzFDvadYs2GKL2BEJqEUgkrL27UMSmDUr\nLD6bPz92RPlj6VK48MKwf/RFF8H48UoCuUSJQKQKNtww1L/v3Bk6dgyzi2T13EN3Wps28O23MGdO\nqOmksYDconUEItU0cWKoWnriiaGraJ11YkeUW956K7QCFiyAAQNg331jR1T4tI5AJMsOPjjMef/w\nQ9h9d3jzzdgR5Yb//Q/OPz9UdT3qKHjtNSWBXKdEIFIDzZvDE09A797hD98114QVssVoxQq48cZQ\nxK+sDObODTvAaVA99ykRiNSQGfToEb75zpkDbduGGjnForQUBg0KZaJfeQVefBHuuSckSckPGiMQ\nSbMxY0LXyF57wQ03hGJ2hai0FIYODWsCtt0WrrsOOnWKHVVx0xiBSI7o1i2smt1xx7C37uWXw5Il\nsaNKn+++g7vuCpv6PPIIDB8OkycrCeQzJQKRDGjYMIwXvPkmfPVV+MZ87bVhy8V89fnn0LdvaOFM\nnQojR4YEsPfesSOTmlIiEMmgzTYL9Ypmzgx/SLfbDq64Aj79NHZkqSkrg0mToHv3UJZ78WJ46aWw\nNkAtgMKhMQKRLPrkk1C36KGHwpTKXr3CbKNaOfaVbO5cGDECHn4YGjeGc84J6yUaNYodmaxJdccI\nlAhEIvjuu/CHduDAMO/+uOPghBPCJu0xVt26h1lPY8eGjWEWLQrxnHxymAWllcD5Ie8SgZnVBl4B\n5rn7ERWeUyKQovHWW6G/fdQo+OmnsFDt4IPD/smNG2fuvJ98Evr6p04Nff3rrgtHHAFdu4bKoLVr\nZ+7ckhn5mAj6AB2BRu7etcJzSgRpVFJSQkJbPqVNpj5P99AlM3Fi+HnxxVCYbbfdwuyjHXaAbbYJ\nJZtTXaTlHmYsffxxWAH9+uvhm//s2WFTmM6dw25g++8P22+f/W/+ujbTq7qJIMqaPzNrARwK/B/Q\nJ0YMxUS/bOmVqc/TLBRna9MG+vQJ8/TnzAmLtF55BZ58MvwxX7gw7JzWtCk0axZaDeXHGL77Lvzx\nX7Ik1PkpKwszfbbeGnbZJZSBbts2PBa7y0fXZm6Itfj7NuBSYP1I5xfJeXXrQrt24adnz18e/+GH\nsHXm4sXhZ8mS8M0fwr/rrQdNmoSfjTYKCSP2H3zJbVlPBGZ2OLDQ3WebWSLb5xfJd+usAy1bhh+R\ndMj6GIGZ9QdOAVYC6xJaBU+4+6nljtEAgYhINeTVYDGAmXUGLqk4a0hERLInF5ax6Nu/iEhEObmg\nTEREsicXWgSY2TFmNsfMfjKzDms4rouZvWNm75vZ5dmMMZ+YWTMzm2xm75nZJDNrsprjPjGzN8xs\ntpm9nO04c1kq15qZ3Zl8/nUza5/tGPPJ2j5PM0uY2TfJa3G2mfWNEWc+MLPBZrbAzFa7J15Vr82c\nSATAm8BRwAurOyC5EvluoAuwI3CCme2QnfDyzhXAZHdvBTyXvF8ZBxLu3t7dd89adDkulWvNzA4F\ntnX37YA/AQOzHmieqMLv7tTktdje3a/PapD5ZQjhs6xUda7NnEgE7v6Ou7+3lsN2Bz5w90/cvRQY\nBXTLfHR5qSswNHl7KHDkGo7VDPPfSuVaW/UZu/tMoImZbZzdMPNGqr+7uhZT4O7TgMVrOKTK12ZO\nJIIU/Q74rNz9ecnH5Lc2dvcFydsLgNVdBA48a2avmNlZ2QktL6RyrVV2TIsMx5WvUvk8Hdgr2ZUx\n3sx2zFp0hafK12bWFpSZ2WRgk0qeusrdx6XwFhrVLmcNn+fV5e+4u69hXcbv3f1LM9sQmGxm7yS/\nbRS7VK+1it9gdY1WLpXPZRawubsvN7NDgNFAq8yGVdCqdG1mLRG4+4E1fIvPgc3L3d+ckOmK0po+\nz+RA0ibuPt/MNgUWruY9vkz+u8jMniI04ZUIUrvWKh7TIvmY/NZaP093/7bc7X+Z2QAza+buebyn\nWzRVvjbnJjRhAAADfElEQVRzsWtodf2ErwDbmVlLM6sHHAeMzV5YeWUscFry9mmEb1e/YmYNzKxR\n8nZD4CDCoL2kdq2NBU4FMLNOwJJy3XHya2v9PM1sY7NQEcnMdidMbVcSqJ4qX5uxis79ipkdBdwJ\nNAeeMbPZ7n6ImW0GDHL3w9x9pZmdB0wEagMPuPvbEcPOZTcCj5rZmcAnwLEA5T9PQrfSk8nfvTrA\nCHefFCfc3LK6a83M/px8/p/uPt7MDjWzD4BlQI+IIee0VD5P4I/AOWa2ElgOHB8t4BxnZiOBzkBz\nM/sMuBaoC9W/NrWgTESkyOVi15CIiGSREoGISJFTIhARKXJKBCIiRU6JQESkyCkRiIgUOSUCEZEi\np0QgIlLklAikYJjZFDM7qMJjvc1swBpe810W4rrAzOaa2fBMn0ukOpQIpJCM5LelCY4DHl7Da7Kx\ntP4c4A/ufkoWziVSZUoEUkieAA4zszoAZtYS2Mzdp5vZ6OS+C29VtvdCsiDam+XuX2Jm15a7f7KZ\nzUxuo3ivmf3md8fM+pjZm8mfC5OP3QtsDUwws97p/g+LpENOFJ0TSQd3/zq59/KhhAqMxwOPJJ/u\n4e6Lzaw+8LKZPe7ua9rlaVVLIbmt4rHAXu7+U7Kr6SRgeLljOgKnE0p51wJmmlmJu59tZgcTtgT9\nTTXNZO395oSywU8By939v9X8CESqRS0CKTTlu4eOS94HuNDMXgNeIvzR3a4K73kA0BF4xcxmA/sD\nW1U4Zm/gSXdf4e7LgCeBfdf0pma2PXCauw8H7gWuAjpUIS6RtFCLQArNWOA2M2sPNHD32WaWIPwx\n7+Tu35vZ88C6FV63kl9/Mapf4fmh7n7VGs7r/HovDWPt4w+nASNgVWtmN2DQWl4jknZqEUhBcffv\ngOeBIfwySLw+sDiZBFoDnSp56QJgIzNrZmbrAIfzyx/y54A/Jrf0JHnMFhVePw040szqJzf6OZK1\n7/ZWD/g0+Z4NgGXu/kIV/rsiaaEWgRSikYSumWOT9ycAZ5vZXOBdQvfQzxzA3UvN7G/Ay4Rt/eau\nOiBsotIXmJQcJC4FziX5Rzx5zGwzezD5eggbAL1e/hyVGAR0NbPNk8fMMLM/uvvj1ftvi1SPNqYR\nESly6hoSESlySgQiIkVOiUBEpMgpEYiIFDklAhGRIqdEICJS5JQIRESKnBKBiEiR+38i3uDq/XoA\n+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1088a5050>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Use exact algebra to find the result for comparison\n",
    "import numpy.linalg as la\n",
    "evals, evecs = la.eigh(Hmat3)\n",
    "print \"Ground state energy: \",evals[0]\n",
    "print \"Ground state wavevector: \",evecs[:,0]\n",
    "\n",
    "# Now set up the simple parameter scan\n",
    "n_alpha = 101\n",
    "alpha_vals = np.linspace(-1,1,n_alpha)\n",
    "energy3 = np.zeros(n_alpha)\n",
    "i=0\n",
    "e_min = 1e30\n",
    "alph_min = 0.0\n",
    "for alpha in alpha_vals:\n",
    "    psi = basis_array[0] + alpha*basis_array[1]\n",
    "    H_psi = -0.5*(d2basis_array[0] + alpha*d2basis_array[1])\n",
    "    add_pot_on_basis(H_psi,Vdiag2,psi)\n",
    "    norm = integrate_functions(psi,psi,num_x_points,dx)\n",
    "    #print alpha, norm\n",
    "    #print Hmat2[0,0] + Hmat2[1,0]*alpha + Hmat2[0,1]*alpha + Hmat2[1,1]*alpha*alpha\n",
    "    energy3[i] = integrate_functions(psi,H_psi,num_x_points,dx)/norm\n",
    "    if energy3[i]<e_min:\n",
    "        e_min = energy3[i]\n",
    "        alpha_min = alpha\n",
    "    i=i+1\n",
    "print \"Min E and alph: \",e_min, alpha_min\n",
    "pl.plot(alpha_vals,energy3)\n",
    "pl.xlabel(r\"Value of $\\alpha$\")\n",
    "pl.ylabel(\"Energy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that this rather simple variational approach does a remarkably good job: we get within 1% of the ground state energy.  Looking at the coefficients of the exact solution, we can see that we'd really need to include further states to fully map out the effects of this rather strong field.  There are various approaches to multi-parameter minimisation which I will describe in a separate notebook."
   ]
  },
  {
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}