TL;DR: The external electric field through the dielectric causes small displacements on the bound charges of the dielectric. These accumulate as an induced surface charge on the dielectric/vacuum interface. The total electric potential is then given as a sum of the external potential, the induced potential of the induced point-like dipole, and the induced potential of the induced surface charge at the dielectric/vacuum interface as given by the multipole rules governing the Poisson equation.
The longer explanation:
I think I can answer the question where the forms of these potentials come from, but be slightly cautious (especially on nomenclature) as by background is in Density Functional Theory. We solve electrostatic equations like this constantly (Hartree-potential) in a polarizable medium (electrons), but in our world, everything is microscopic and nicely continuous.
As the example PDF assumes, there are no charges inside the cavity $0 < r < R$, where the Eq. (1) is fulfilled. And it also is true that there are no free charges outside the cavity, r>R and the Eq. (1) is also fulfilled. What particular detail it omits (which might lead to confusion) is that there will be induced charge shell at $r=R$. The differential equation can be solved with mathematical constraints, but I think it is crucial for the physical picture to understand the induced charges.
When there are extra charges in the system, the potential can be solved using following equation, so called Poisson equation
$$ \nabla^2 \phi = n,$$
where n is the total charge density. This can be divided into external, and induced charge densities
$$ \nabla^2 \phi = n_{ext} + n_{ind}. $$
The induced charges are due to the bound charges in the dielectric medium polarizing (i.e. moving away from their neutralizing background charge), or due to the dipole moment of the point-like object. In the case of the surface charge, the width of this charge will be infinitesimal, but the charge density on this charged shell diverges too to keep everything in control. First, let's verify that this is indeed the case. Let's take the left hand side of the Poisson equation (the Laplace operator) and plug in the already solved potential. I will omit LaTeX, and do this in Maple, as it gets rather tedious.
First, let's define the potential they solved
> F:=piecewise(r<R, p*cos(theta)/r^2 - 1/(1+2*eps)*(2*(eps-1)*p/R^3+3*eps*Einf)*r*cos(theta), r>=R, 1/(1+2*eps)*(3*p-(eps-1)*R^3*Einf)*cos(theta)/r^2-Einf*r*cos(theta));
{ /2 (eps - 1) p \
{ |------------- + 3 eps Einf| r cos(theta)
{ | 3 |
{ p cos(theta) \ R /
{ ------------ - ----------------------------------------- r < R
{ 2 1 + 2 eps
{ r
{
{ 3
{ (3 p - (eps - 1) R Einf) cos(theta)
{ ------------------------------------ - Einf r cos(theta) R <= r
{ 2
{ (1 + 2 eps) r
We can verify that the potential is continuous at R.
> simplify(limit(F,r=R,right)-limit(F,r=R,left));
0
Now we can apply the Laplacian operator (in spherical coordinates) to the potential.
> simplify(convert(1/r*diff(r*F,r$2)+1/(r^2*sin(theta))*diff(sin(theta)*diff(phi,theta), theta), piecewise,r));
{ 0 r~ < R~
{
{ undefined r~ = R~
{
{ 0 R~ < r~
We indeed find, that the induced density is zero in the regions described by the example, and Eq. (1) can be used. However, Maple cannot immediately do the $r=R$ case. This is because already the gradient of this potential is discontinuous. Let's see what the discontinuity in the gradient looks like. There is no $\phi$ dependence, so gradient discontinuity to $\phi$ direction will be 0. There is $\theta$ dependence, but we get:
> limit(diff(phi,theta)/r,r=R,right)-limit(diff(F,theta)/r,r=R,left);
0
which was one of the boundary conditions given by the example, and for gradient discontinuity in $r$ direction we get
> simplify(limit(diff(F,r),r=R,right)-limit(diff(F,r),r=R,left));
3
3 cos(theta~) (Einf~ R~ + 2 p~) (eps~ - 1)
-------------------------------------------
3
(1 + 2 eps~) R~
Note, that this also has the $\cos(\theta)$ symmetry. Indeed, there is a discontinuity in the gradient of the potential i.e. the electric field. Discontinuity in the radial component of the electric field means that there is surface charge at $r=R$.
Now, that we understand the form of induced charge in the system, understanding the forms of chosen potentials is easier. If you look at this link, you will see how you can write the electrostatic potential as an integral over the Coulomb kernel $\frac{1}{|r-r'|}$. There exists a particular expansion for this potential in spherical symmetry, given with r< and r> terms, where $r\leq \min(r,r')$ and $r\geq \max(r,r')$. $Y_lm$ are the spherical harmonic functions, and all we need to know that $A\cdot\cos(\theta)$ is one of the $l=1, m=0$ spherical harmonic, and these functions are orthogonal. Therefore the angular integrals all vanish except the $l=1,m=0$ term, where it yields a constant (omitted, and so is the angular part of the integral for brevity).
Now, finally we can get to the forms of these potentials. First, let's look at the case r>R and r'<=R. This means, let's look at the potential outside the sphere, induced by charge inside the sphere (and the surface).
We can write
$$ \phi_{out,due~in}(r) = \frac1{r^2} \cos(\theta) \int_{0}^{R+} dr' r' n_{ind}(r') $$.
The integral integrates the 'dipole moment' of the density inside the sphere. And since we multiply this integrated constant by 1/r^2 cos(theta), it shows outside the sphere exactly the same way as a point like dipole. The integral of dipole moment, $p_{out}$ will be sum of the point dipole dipole moment at the center, $p$, plus the induced dipole moment on the dielectric surface $p_{ind}$:
$$ \int_{0}^{R+} dr' r' n_{ind}(r') = p_{out} = p + p_{ind}. $$
Let's then analyse the case where r>R, r'>R, but r' is very large. We can write
$$ \phi_{out,due~far~away} = r \cos(\theta) \int_{large}^{\infty} dr' n(r') / r'^2. $$
We see that the linear electric field in the problem, E_inf, could be pictured to be cause by an external charge density with $cos(\theta)$ symmetry far away, integrating to $E_{\infty} = \int_{large}^{\infty} dr' n(r') / r'^2$.
Let's then analyse the case where $r'=0$, and $0<r<R$. We have the point like dipole at $r'=0$, which can be described if we allow distribution like induced density.
$$ \phi_{in,due~dipole} = \frac{1}{r^2} \cos(\theta) p $$,
thus we get what we expect. The dipole moment integral over a dipole moment yields naturally the dipole moment.
And then finally, let's analyse the case where $0<r<R$, and $r'>R$.
$$ \phi_{in,due~out} = r \cos(\theta) \int_{R}^{\infty} dr' n_{ind} / r^2 $$.
Here the integral $\int dr' n_{ind} / r^2$, gives $E_{in}$.
The Poisson equation is additive, so we can sum the in due dipole and in due out, to get
$$ \phi_{in} = \frac{1}{r^2} \cos(\theta) p + r \cos(\theta) E_{in} $$
and
$$ \phi_{out} = r \cos(\theta) E_{\infty} + \frac{1}{r^2} \cos(\theta) (p + p_{ind}) $$.
This concludes, why the potentials in and out, look the way they look.
What we actually have come up against, is the powerful theory of multipoles in a Laplacian equation. Given a Laplace equation, we notice, that if we know certain integrals of charge density shells with domain outside of our own shell, and certain integrals of the inside spherical shells of our location (called multipole moments), we can solve the potential at our location.
And one more remainder, that the integrals I present are only sketched to give the idea, they miss factors and part of angular integrals etc.