By this point in the sequence of tutorials, we have developed enough skills with Mathematica™ to construct a relatively-full example that is a classic in quantum mechanics: the harmonic oscillator. In the development below, we will be assuming that we have already derived the wavefunctions and energy levels for this system (we can certainly derive them in Mathematica™, and I may add a tutorial on how to do that at a later point). Over the course of this tutorial, we will generate the figure shown to the right, showing the probability density functions and energy levels for the harmonic oscillator for levels v=0 to 10.

Let’s start with the potential energy curve. We are also going to set up some specific conditions that we will use for making plots: the force constant will equal 1, the mass will equal 1, and reduced Planck’s constant will equal 1 (*i.e.*, we are operating in atomic units). Notice the usage of Set, the definition of a function with a variable argument (in this case, the x coordinate), the replacement rules that allow for the application of conditions temporarily, a two-dimensional plot of a function implementing those conditions (with a bunch of aesthetic modifiers that you don’t need to worry about, but which will make our final result easier to visualize), and the fact that we can assign the plot itself to a variable for later reference.

Now let’s define the wavefunctions. There are a couple of advanced techniques being used in that definition below, which I will explain, but hopefully you will be able to see how you can do the same thing without using those advanced techniques. First, let me point out the aspects of the code that we went over in earlier tutorials:

- The use of Set to define the function. Since this particular function is not well-defined when the quantum number v is not an integer, you could argue that one should define this function using SetDelayed instead. The use of ConditionalExpression (described below), however, means that the function is only defined in cases where v is in the correct domain.
- The use of Greek letters to make the code follow standard chemistry notation more closely.
- The separation of one function parameter (the quantum number) as a subscript from the other function parameter (the x-coordinate) as being within the brackets.
- The use of replacement rules to apply our specific conditions to the wavefunction (in the last line).

The advanced techniques that you can do without if you so choose are:

- ConditionalExpression. This says that the part before the comma is the value of the function, but only under the conditions that are listed after the comma. As I mentioned above, you can leave this out and use DelayedSet instead.
- Module. This is a way to define “local” variables for a set of statements. Here, α and y are only being defined within the Module. We assign values to them, and then use those values in the ConditionalExpresson to simplify the form of the function. Outside of the Module, those two local variables are not modified. This is a convenient way to ensure that you don’t litter the kernel with definitions that might interfere with later work you are doing.

We also need the energies. Notice, in the code below, the use of a Table to print out some test values to make sure the expression works, and the use of replacement rules to substitute in our conditions.

What we would really like to do is to plot the energy levels on the potential energy curve. We have the y values (the energies), but we also need to know the x-values to limit the horizontal span of the lines. The x values are the places where the energy level are equal to the potential energy. Let’s see how to do that.

It is instructive to consider the first line above, and how it arose as I was coding this. Many students look at that kind of line and think, “Wow, there’s a lot going on at once… how do you put it all together right away?” And the answer is that you don’t. You do it stepwise, making sure that each part works first. Let’s walk through how I did it.

First, I thought, “The basic thing I want to do is find the x values where the energies are equal. That sounds like a Solve statement. Let’s pick v=3 just to have something to work with:

Ok, that’s looking like a good start. But let’s simplify our life by substituting in for the constants. Copy-paste and modify:

Notice the use of parentheses to ensure that the replacement rule “/.conditions” applies to both sides of the equation. Ok, this is good so far, but I don’t really want the result to ** be** a replacement rule, so we should

**the replacement rule:**

*use*That says, “give me x, where x has been replaced with these two values.” So let’s put that in conjunction with our earlier statement.

Wonderful. Now let’s generalize this so that we can do it at any quantum number, and we will write this as a function where the only argument is a subscript:

And we’re done. Now in reality, I did not actually keep all of those lines. I started with the first one, checked it, and then made sequential modifications to it until it worked like I wanted. Like this:

Ok, so we now have the crossings. We want lines running from {leftCrossingx,energy} to {rightCrossingx,energy}. Let’s do see what we can do. In particular notice the use of the Table within the Graphics statement. As with many other multi-component lines, this one was built up by working from the inside out, similar to my extended example above.

Now it would be nice to plot the wavefunctions on each line. Or actually, the wavefunction squared would probably be better. To do that, we need to offset each wavefunction vertically from zero by the magnitude of the energy level. And just from a graphics prettiness standpoint, we should scale them slightly. Fortunately, that is straightforward.

Now we just combine this with the energy level diagram, and we are done.

Our final result is a professional-looking plot that contains a huge amount of information.