Solving ODE using python

Please find the mistake in this code

Thanks!

Question 12: Consider the one-dimensional, time-independent Schroedinger equation describing a particle with wave function (x) and energy E : 2m2dx2d2(x)+V(x)(x)=E(x). Values for the constants m, and others required for the calculation in SI units are given in the code cell below. (Note: the energy unit eV is taken as elementary charge e times Volts with e given as a constant.] Tasks: (a) Consider a square potential well with V(x)=0 in between two infinitely high walls separated by a distance L equal to the Bohr radius, i.e. for all x in the interval [0, L ]. Write a function solve(energy, func) which takes the parameter energy and a python function f unc. This function should solve the Schroedinger ODE for the case described above and return only the final value of (L) at the boundary L. Write a function eqn ( x,y, energy) containing the two first-order differential equations equivalent to the Schroedinger equation. Use this function as input for the solve(energy, func) function. This detailed naming scheme is required for the first tests of your code. charge value). For the initial condition, see technical hint below. For the number of values to solve for (x) the value 1000 is recommended. The required precis the reas finding method is 1.e22 (set explicitly since less precision has consequences for part (c)). The result of your calculation should be a number between 134eV and 135eV. One of the tests will assess your solve(energy, func) function for a distorted potential well. (b) Consider the harmonic potential V(x)=V0x2/a2, where V0 and a=1011m are constants. Limit the infinite range of the variable x to the interval [10a,10a] with V0=50eV. The harmonic oscillator is known to have equidistan energy eigenvalues. Check that this is true, to the precision of your calculation, by calculating the ground state and the first 2 excited states. (Hint: the ground state has energy in the range 100 to 200eV. You may wish to program a different solve2(energy) function to part (a) for this exercise.) In order to test your result, write a function result() which must return the difference of calculated energy eigenvalues in eV as a single number. As stated above, the absolute difference between the ground and first state and the first and second state should be very close to equal. Either of those two numbers should be returned. Consider that this is merely a function to enable automated testing and hence should really only return the requested value. Calculations should take place earlier in your code so that result() can run without an input argument. (c) Plot the normalized wavefunctions, for the ground state energy and the two excited state energies as obtained without normalization in part (b). Note that normalization of a wavefunction requires the calculation of the integral of the square of the wavefunction (x)2dx and subsequent division of (x) by the square root the integral. Provide a plot title and appropriate axis labels, plot all three wavefunctions onto a single canvas using color-coded lines and also provide suitable axis limits in x and y to render all curves well visible. Technical Hint: This is not an initial value problem for the Schroedinger ODE but a boundary value problem! This requires additional effort, a different method to the previous ODE exercises and is hence an 'unseen' problem. Take a simple initial condition for both problems at x0=0 or x0=10a, respectively: (x0)=0 and d(x0)/dx=1. Use that as a starting point to solve the ODE and vary the since the first boundary condition is already satisfied due to the choice of initial condition. This requires an initial guess for the energy interval in which the solution might be and a computational method for root finding. Search scipy for root finding methods and select one which does not require knowledge of the derivative. Root finding is appropriate here since the boundary condition to satisfy is (x)=0 to a given precision. Another hint would be to consider the required precision since in SI units, absolute values for E are very small indeed. Technical Hint 2: part (a) requires using solve_ivp in order to solve the boundary value problem in combination with a root finder function. The sciPy boundary solve_bvp , can be used for parts (b) and (c), if you wish. The root finding, solve_ivp method will also work, however. Quantum physics background For those lacking the first year quantum module, the boundary condition for both exercises is that (x)=0 at each potential boundary. These are only fulfilled at specific, discrete energy values E, called energy eigenvalues, where the smallest of these is called the ground state energy. \# axis limits plt. xlim(10L,10L) plt.ylim( (2,2) \# legend plt. legend() plt.show()