For this example lets find the wavefunction of the bound states of a quantum particle in a well using scaled units where:

First lets solve the Schrodinger equation with `DSolve`

:

```
(* Region 1 and 3 *)
DSolve[psi''[x] - d^2 psi[x] == 0, psi[x], x]
(* Region 2 *)
DSolve[psi''[x] + k^2 psi[x] == 0, psi[x], x]
```

Two of these solutions are not physical as the wavefunction must be normalizable, and one of the integration constants can be absorbed into the unknown value for the energy. The wavefunctions, and their derivatives, in the three regions are then:

```
psi1 = Exp[d x]
d1 = D[psi1, x]
psi2 = a2 Cos[ k x] + a3 Sin[ k x]
d2 = D[psi2, x]
psi3 = a4 Exp[ - d x]
d3 = D[psi3, x]
```

In order to do the calculation numerically, we need to set the well depth (20) and width (1) in scaled units.

```
wellDepth = 20;
wellWidth = 1;
```

Next, we set up a system of equations that match the boundary conditions of both the wavefunction and its derivatives at the boundaries. Use a replacement rule (`/.`

) to substitute the energy (en) and the well depth for the corresponding values of k and d. This linear system has 4 equations and 3 unknowns, so the `Eliminate`

function can reduce the system to one equation:

```
eqs = {psi1 == psi2 /.x->0, d1 == d2 /.x->0, psi3 == psi2/.x->wellWidth, d3 == d2/.x->wellWidth}/.{d->Sqrt[-en], k->Sqrt[wellDepth+en]};
MatrixForm[eqs]
eq1 = FullSimplify[Eliminate[eqs, {a2, a3, a4}]]
```

Finding the roots of the equation is made easier by solving for zero and plotting. This shows that the roots are at approximately - 15 and - 4. Then use `FindRoot`

to find the roots of this equation. Note that you will need to supply a guess to `FindRoot`

as a starting place of the search for the root.

```
Plot[ en (20 + en) + 100 Sin[Sqrt[20 + en]]^2, {en, -wellDepth, 0}]
e0 = en /. FindRoot[eq1, {en, -15}]
e1 = en /. FindRoot[eq1, {en, -4}]
```

`FindRoot`

can also find zeros of system of equations by supplying all of the unknowns for the system, as well as a guess for each unknown. Finding the constants for the ground state yields

```
roots0 = FindRoot[ eqs, { en, e0}, {a2, 1}, {a3, 1}, {a4, 1}]
```

Solve for the $0^th$ wave function (wf0) in each of the 3 regions and plot the result:

```
wf0 = ( { psi1, psi2, psi3} /. {d -> Sqrt[-en], k -> Sqrt[wellDepth + en]}) /. roots0
psi0Plot = Plot[ If[x < 0, wf0[[1]], If[x < 1, wf0[[2]], wf0[[3]]]], {x, -1, 2}]
```

Finding the constants for the first excited state yields

```
roots1 = FindRoot[ eqs, { en, e1}, {a2, 1}, {a3, 1}, {a4, 1}]
```

Solve for the $1^st$ wave function (wf1) in each of the 3 regions and plot the result:

```
wf1 = ( { psi1, psi2, psi3} /. {d -> Sqrt[-en], k -> Sqrt[wellDepth + en]}) /. roots1
psi1Plot = Plot[ If[x < 0, wf1[[1]], If[x < 1, wf1[[2]], wf1[[3]]]], {x, -1, 2}, PlotStyle -> Red]
```

Finally put it all together with a picture of the well and the wavefunctions

```
wellPic = Graphics[{Line[{{0, 0}, {1, 0}}], Line[{{0, 0}, {0, 2}}],
Line[{{1, 0}, {1, 2}}], Line[{{-1, 2}, {0, 2}}],
Line[{{1, 2}, {2, 2}}], Text["\[Psi]1", {-0.3, 1}],
Text["\[Psi]2", {0.5, 1}], Text["\[Psi]3", {1.3, 1}]}];
Show[wellPic, psi0Plot, psi1Plot]
```