Question: How do I get Maple to fsolve using the Shooting method for boundary value problems?

This worksheet covers the energy levels of a muon bound to a lead atom (Pb-206). All plots have been commented out and inserted as pictures. First define the variables for this atom, and the muonic mass ms in units of the free electron mass:

where Z is the atomic number and A is the mass number. The Z and A values can be changed for any other atom or isotope. Modelling the atom as hydrogenic, we have the equation

where n is the orbital occupied by the muon and En is the energy in Hartree units. The average radius a of the ground state orbit in Bohr is

This is much smaller than the Bohr radius due to the muonic mass being bigger by 2 orders of magnitude and the higher nuclear charge. The atomic radius of the lead atom is determined by

This is over twice the above average orbital radius, meaning that the ground state orbit lies within the nucleus. This means we cannot accurately treat this atom as hydrogenic, since the muon will only be attracted towards the portion of the nuclear charge closer to the centre. Instead, we must use a model of nuclear structure that accounts for the muon being inside the nucleus, treating the nucleus as a spherical charge distribution. Solving the radial Schrodinger equation for the sphere using a uniform spherical charge distribution, then applying Gauss's law for the charge enclosed by the nuclear surface gives a potential distribution which can be modelled using a piecewise function for the different potentials inside and outside the nucleus:

While the point charge model unrealistically goes to an infinite potential as the muon gets closer to the centre of the nucleus, the spherical model instead shows it tailing off inside the nucleus as more of the nuclear charge is further from the centre than the muon itself. To model the bounds of the energy of the ground state 1s, we use the Shooting method for boundary value problems.

We are looking for a solution with l=0, and the offset r0 to prvent division by 0 must be much smaller than the Bohr radius given the scale of the orbit. We need to estimate bounds for the energy and set an upper bound for the effective range of the potential. This upper bound must be low enough for Maple to calculate the function, which unfortunately means only a few Bohr radii out.

These bounds show that the ground state binding energy is between -39 and -38 Hartrees, since fnl tends to negative infinity at one bound and positive infinity at the other. To home in on the correct value, we plot the function between the bounds and find the zero.

The plotted function tends to a straight line with no nodes, so is the solution to 1s. This gives a binding energy of -11.94 Hartrees for the ground 1s state. The 2p state is the lowest with l=1, so the same method using l=1 will find the energy of the 2p state. First find upper and lower bounds:

Then solve the differential equation: