This is my solution to problem #3 on the Midterm, Spring 1997.

Statement of the Problem

Consider the expression from Jackson's problem 1.5 for the potential of a neutral hydrogen atom:

      phi = q exp(-alpha*r)/r (1 + alpha*r/2)

Using MatLab and following the general lines of Jim Lockhart's example, explore the charge distribution corresponding to this potential.

Solution to the Problem

Okay, as suggested, I split the equation for phi into two parts. The first part is:

      phi1 = q * exp((-a).*r) ./ r

The second part is:

      phi2 = q * exp((-a).*r) ./ r .* a .* r ./ 2

The code for the matlab program, midterm3.m reveals that q=1, and a=2/a0=2.

To begin, let's look at the laplacian of phi1. The graph reveals a double spike near the origin. It would make sense to have a spike near the origin since one the laplacian of (1/r) yields a delta function. I do admit that it does not seem quite right to have a double spike. I would expect a single spike. I wonder if it has to do with the fact the my matrix is being divided by zero at the origin. A new graph plotted excluding the origin reveals a single spike. The compliation of this new graph showed no warnings from matlab of a division problem. Continueing on, the graph of the laplacian of phi2 shows a spike (not as sharp of a spike as the last graph) but in the opposite direction. This also makes sense since this is the part due to the electron. The other part, phi1, is due to the proton. Therefore, the two parts should be in opposite directions. Onward to the graph of the laplacian of phi. This is the sum of phi1 and phi2. The graph of phi does not look much different from the first graph (phi1). Hmm, I expected this graph to be different. This is because the graph of phi is the the addition of phi2 to phi1. The resultant graph just does not seem to be totally correct. Since I cannot find the discrepancy, I will continue on as though the graph is fine. By the way, if you are interested in looking at the graph from a different perspective, click here.

Since rho, the charge distribution, is: rho = (the laplacian of phi)/(-4*Pi), I simply multiplied the laplacian of phi by the appropriate constants. A graph of rho reveals a spike in the negative direction. A view of the graph from the side is also interesting. Checking the laplacian of phi matrix to see how close numerically it is to satisfying the equation laplacian of phi = -4 pi rho, it seems that they are very close but not exact. How large the discrepancy is depends upon where on the x-y-z grid one looks and to what precision one is interested. The graph seems to become less precise around the origin. Oh, if you are interested in the code to the matlab program, click here.