%% Problem Set J, Problem 3, Parts (a) & (b), but only subpart (ii)

clear all
close all

%
%% (a) Gauss' Law
% Let g(r) be the gravitational field strength, rho(r) the mass density,  gravconst the gravitational constant. Then we have

syms a r z theta h rho(r) gravconst real
g = @(a) 4*pi*gravconst*int(int(int(rho(r)*r, r, 0, a), z, 0, h),...
    theta, 0, 2*pi)/(2*pi*a*h);
g(a)

% We use this expression below.
%
%% (ii) Decrease like the 4th power

gravconst=1; rho = @(r) 1/(1+r^2)^2
gravfieldstrength = simplify(4*pi*gravconst*int(r*rho(r), r, 0, a)/a)
fplot(gravfieldstrength, [0 5])

%
% The field strength increases initially -- apparently linearly for a while -- reaches a maximum at about r = 1, then begins to decay, rather slowly it would appear, although the equations suggest in inverse proportion to the distance. This is rather slow decay, so the field is still substantial far out from the mass. Here's another graph that substantiates the claim.

figure; fplot(gravfieldstrength, [0 15])

%
%% (b) Poisson's Equation
% We use the Laplacian in cylindrical coordinates (actually, it is essentialy polar coordinates here since no z is involved in the computation). In fact \nabla^2(\rho(r)) = \rho''(r) + (1/r)\rho'(r) (also using that the derivative with respect to \theta is also zero).
%% (ii)

syms r real; rho = @(r) 1/(1+r^2)^2
syms t phi(t) real;
sol = dsolve(diff(phi(r),r)/r + diff(phi(r),r,2) == ...
    4*pi*gravconst/(1 + r^2)^2, r);
simplify(sol)

% The constant in front of the \log r must match 2\pi and so the term we must deal with is const + \pi \log(1+ r^2), whose gradient is proportional to r/(1+r^2), which again is what we got in part (a).