Verify $y=x^aZ_p\left(bx^c\right)$ is a solution to $y''+\left(\frac{1-2a}{x}\right)y'+\left[(bcx^{c-1})^2+\frac{a^2-p^2c^2}{x^2}\right]y=0$ In order for the question that I have to make any sense I must first include some background information as given in my textbook:

The standard form of Bessel's differential equation is $$x^2y^{\prime\prime}+xy^{\prime} + (x^2 - p^2)y=0\tag{1}$$ where $(1)$ has a first solution given by $$\fbox{$J_p(x)=\sum_{n=0}^\infty\frac{(-1)^n}{\Gamma(n+1)\Gamma(n+1+p)}\left(\frac{x}{2}\right)^{2n+p}$}\tag{2}$$ and a second solution given by $$\fbox{$J_{-p}(x)=\sum_{n=0}^\infty\frac{(-1)^n}{\Gamma(n+1)\Gamma(n+1-p)}\left(\frac{x}{2}\right)^{2n-p}$}\tag{3}$$ where $J_p(x)$ is called the Bessel function of the first kind of order $p$. 
Although $J_{−p}(x)$ is a satisfactory second solution when $p$ is not an integer, it is customary to use a linear combination of $J_p(x)$ and $J_{−p}(x)$ as the second solution. Any combination of $J_p(x)$ and $J_{−p}(x)$ is a satisfactory second solution of Bessel’s equation. The combination which is used is called the Neumann (or the Weber) function and is denoted by $N_p(x)$ where $$N_p(x)=\frac{\cos(\pi p)J_p(x)-J_{-p}(x)}{\sin(\pi p)}\tag{4}$$
Full details on the derivation of $(2)$ as a solution to $(1)$ can be found here in my previous question.
Many differential equations occur in practice that are not of the standard form $(1)$ but whose solutions can be written in terms of Bessel functions. It can be shown that the differential equation: $$\fbox{$y^{\prime\prime}+\left(\frac{1-2a}{x}\right)y^{\prime}+\left[\left(bcx^{c-1}\right)^2+\frac{a^2-p^2c^2}{x^2}\right]y=0$}\tag{5}$$ has the solution $$\fbox{$y=x^aZ_p\left(bx^c\right)$}\tag{6}$$ where $Z_p$ stands for $J_p$ or $N_p$ or any linear combination of them, and $a,b,c,p$ are
  constants.
To see how to use this, let us “solve” the differential equation: $$y^{\prime\prime}+9xy=0\tag{7}$$ If $(7)$ is of the type $(5)$, then we must have $$1-2a=0$$ $$2(c-1)=1$$ $$(bc)^2=9$$  $$a^2-p^2c^2=0$$ from these $4$ equations we find
  $$a=\dfrac12$$ $$c=\dfrac32$$ $$b=2$$ $$p=\dfrac{a}{c}=\dfrac13$$
Then the solution of $(7)$ is $$y=x^{1/2}Z_{1/3}\left(2x^{3/2}\right)\tag{8}$$ This means that the general solution of $(7)$ is $$y=x^{1/2}\left[AJ_{1/3}\left(2x^{3/2}\right)+BN_{1/3}\left(2x^{3/2}\right)\right]\tag{9}$$ where $A$ and $B$ are arbitrary constants.


Finally,$\color{#180}{\text{ my goal is to show that }}$${(6)}$ $\color{#180}{\text{is a solution to }}$$(5)$. 
However to gain some insight, I must first be able to show that $(8)$ or $(9)$ is a solution to $(7)$.

So my attempt goes as follows:
So I need to compute $y^{\prime\prime}$ or at least to begin with, $y^{\prime}$; It is at this point where I am immediately stuck as I do not understand how to differentiate $$y=x^{1/2}Z_{1/3}\left(2x^{3/2}\right)\tag{8}$$ as I'm confused as to how to take the derivative of the $Z_{1/3}\left(2x^{3/2}\right)$ factor.
Could someone please provide some hints or advice on how I would go about carrying out this differentiation?
 A: $$\fbox{$y^{\prime\prime}+\left(\frac{1-2a}{x}\right)y^{\prime}+\left[\left(bcx^{c-1}\right)^2+\frac{a^2-p^2c^2}{x^2}\right]y=0$}\tag{5}$$
Let $y(x)=x^\alpha Y(X)\quad$ with $X=\beta x^\gamma$
$\qquad \alpha, \beta , \gamma$ are constants.
$\frac{dX}{dx}=\gamma\beta x^{\gamma-1} \quad;\quad \frac{d^2 X}{dx^2}=\gamma(\gamma-1)\beta x^{\gamma-2}$
$\frac{dy}{dx}=\alpha x^{\alpha-1}Y + x^\alpha \frac{dY}{dX}\frac{dX}{dx} 
= \alpha x^{\alpha-1}Y + \gamma\beta x^{\alpha +\gamma-1}\frac{dY}{dX}$
$\frac{d^2y}{dx^2}=\alpha(\alpha-1) x^{\alpha-2}Y + 2\alpha x^{\alpha-1}\frac{dY}{dX}\frac{dX}{dx}+x^\alpha \frac{d^2 Y}{dX^2}(\frac{dX}{dx})^2 +x^\alpha \frac{dY}{dX}\frac{d^2 X}{dx^2}$
$\frac{d^2y}{dx^2}=\alpha(\alpha-1) x^{\alpha-2}Y 
+ (2\alpha+\gamma-1) \gamma\beta x^{\alpha+\gamma-2}\frac{dY}{dX}
+ \gamma^2\beta^2 x^{\alpha+2\gamma-2} \frac{d^2 Y}{dX^2} $ 
Putting them into (5) :
$$y^{\prime\prime}+\left(\frac{1-2a}{x}\right)y^{\prime}+\left[\left(bcx^{c-1}\right)^2+\frac{a^2-p^2c^2}{x^2}\right]y = 0 =\\
\alpha(\alpha-1) x^{\alpha-2}Y 
+ (2\alpha+\gamma-1) \gamma\beta x^{\alpha+\gamma-2}\frac{dY}{dX}
+ \gamma^2\beta^2 x^{\alpha+2\gamma-2} \frac{d^2 Y}{dX^2}
+(1-2a)\left(\alpha x^{\alpha-2}Y + \gamma\beta x^{\alpha+\gamma-2}\frac{dY}{dX}\right)
+\left[\left(bcx^{c-1}\right)^2+\frac{a^2-p^2c^2}{x^2}\right]x^\alpha Y$$
With $X=\beta x^\gamma$ and after simplification :
$$ X^2\frac{d^2 Y}{dX^2} 
+ \frac{(2\alpha+\gamma-2a)}{ \gamma} X\frac{dY}{dX} 
+\left[\frac{ b^2c^2 }{\gamma^2\beta^{2c/\gamma}}X^{2c/\gamma} + \frac{a^2-p^2c^2+\alpha(\alpha-2a)}{\gamma^2} \right]Y =0$$
In order to obtain the standard form of Bessel equation : $X^2Z_p''+XZ_p'+(X^2-p^2)Z_p=0$ , the relationships are :
$$\begin{cases} 
 \frac{(2\alpha+\gamma-2a)}{ \gamma}=1 \\
2c/\gamma=2 \\
\frac{ b^2c^2 }{\gamma^2\beta^{2c/\gamma}}=1\\
\frac{a^2-p^2c^2+\alpha(\alpha-2a)}{\gamma^2}=-p^2 \\
Y(X)=Z_p(X)
\end{cases} 
\quad\to\quad 
\begin{cases}
\alpha=a \\
\beta=b \\
\gamma=c
\end{cases}
$$
Thus $y(x)=x^a Y(X)=x^a Z_p(X)\quad$ with $X=b x^c$
Hense, the basic solutions of eq.(5) are $y(x)=x^a Z_p(b x^c)$
A: This answer is too long for a comment and refers to BLAZE's comment above as to how to compute $$\frac{d}{dx}\frac{dY}{dX}$$
Setting the Scene
$y(x)=x^{\alpha}Y(X),\quad X=\beta x^{\gamma}$
Now it is obvious that
\begin{align}
\frac{dy}{dx} &= \alpha x^{\alpha - 1}Y+x^{\alpha}\frac{dY}{dx} \\
              &= \alpha x^{\alpha-1}Y+x^{\alpha}\frac{dY}{dX}\frac{dX}{dx}
\end{align}
We could expand as JJacquelin does but let's leave the calculation here. Now computing the second derivative
\begin{align}
\frac{d^{2}y}{dx^{2}} &= \frac{d}{dx}\left(\alpha x^{\alpha-1}Y+x^{\alpha}\frac{dY}{dX}\frac{dX}{dx} \right) \\
&= \alpha (\alpha - 1)x^{\alpha-1}Y+\alpha x^{\alpha-1}\frac{dY}{dx}+\alpha x^{\alpha-1}\frac{dY}{dX}\frac{dX}{dx} + x^{\alpha}\frac{d}{dx}\left(\frac{dY}{dX}\frac{dX}{dx}\right) \\
&= \alpha (\alpha - 1)x^{\alpha-1}Y+\alpha x^{\alpha-1}\frac{dY}{dX}\frac{dX}{dx}+\alpha x^{\alpha-1}\frac{dY}{dX}\frac{dX}{dx}+x^{\alpha}\frac{d^{2}Y}{dX^{2}}\left(\frac{dX}{dx}\right)^{2}+x^{\alpha}\frac{dY}{dX}\frac{d^{2}X}{dx^{2}}\\
&=\alpha (\alpha - 1)x^{\alpha-1}Y+2\alpha x^{\alpha-1}\frac{dY}{dX}\frac{dX}{dx}+x^{\alpha}\frac{d^{2}Y}{dX^{2}}\left(\frac{dX}{dx}\right)^{2}+x^{\alpha}\frac{dY}{dX}\frac{d^{2}X}{dx^{2}}
\end{align}
I hope this clears things up a little for BLAZE.
Apologies to the moderators with regards to this not being a full answer to the question, but I thought the level of computation was too long for a comment.
A: I will here show how one could go about to derive the solution of your ODE without knowing a priori what the answer should be.

Your ODE $(5)$ can after multiplication by $x^2$ be written
$$x^2y^{\prime\prime}(x)+x(1-2a)y^{\prime}(x)+\left[\left(bcx^{c}\right)^2+a^2-p^2c^2\right]y(x)=0$$
This looks very close to the Bessel ODE so we will try to see if we can bring it closer to that form. The term with $x^c$ on the right-hand side is troublesome so lets try to simplify this bringing it closer to the $x^2$ form we have in the Bessel ODE. This motivates a change of variables $x\to x^c$. Another reason for why this is a natural thing to try is the fact that the derivatives are all on the (logarithmic) form $x^ny^{(n)}(x)$. This means that a change of variables $x\to x^n$ will not change the basic structure of the derivative-terms in the ODE - just the coefficient in front of them. 
We therefore take $z = x^c$ and compute
$$x\frac{d}{dx} = x\frac{dz}{dx}\frac{d}{dz} = cz\frac{d}{dz}$$
$$x^2\frac{d^2}{dx^2} = c^2 z^2\frac{d^2}{dz^2} + c(c-1) z \frac{d}{dz}$$
to get that the ODE transforms into
$$z^2 y''(z) + z\left(1-\frac{2a}{c}\right)y'(z) + [(bz)^2 + \frac{a^2}{c^2} - p^2]y(z) = 0$$
where a prime is now a differential with respect to $z$. The $z^2$ term is now close to the desired form and we just need a simple scaling of the variables $w = bz$ to get it right
$$w^2 y''(w) + w\left(1-\frac{2a}{c}\right)y'(w) + [w^2 + \frac{a^2}{c^2} - p^2]y(w) = 0$$
Now this is starting to look more and more like the Bessel ODE, however we still have the troublesome $-\frac{2a}{c}$ term. There is no simple change of variables that can solve this, but a transformation $y(w)  = w^n g(w)$ could potentially work as this has the effect of modifying the relative coefficents of the $w^nf^{(n)}(w)$ terms. Another way that could have lead us to try this is the fact that if we did not have the $w^2y(w)$ term in the ODE then $y(w) = w^{\frac{a}{c} \pm p}$ would be a solution to the ODE so it would be natural to try to factor out the power-law growth by trying an ansatz on the form $y(w) = w^n g(w)$.
Substituting $y(w) = w^n g(w)$ into the ODE above leads to
$$w^2g''(w) + \left(1-\frac{2a}{c} + 2n\right)wg'(x) + \left[w^2 + \frac{a^2}{c^2} -2n\frac{a}{c} + n^2 - p^2\right]g(w) = 0$$
We see that the choice $n = \frac{a}{b}$ does the job and brings the term $wg'(w)$ to the desired form (relative to $w^2g''(w)$) and as a bonus the last terms also happens to simplify for this choice of $n$ leaving us with
$$w^2g''(w) + wg'(w) + \left[w^2 - p^2\right]g(w) = 0$$
which is exactly the ODE for the Bessel-functions. Backtracking our steps we get the solution in terms of the original variables
$$f(x) = w^nZ_p(w) = (bz)^nZ_p(bz) = b^{\frac{a}{c}} x^{\frac{a}{c}\cdot c} Z_p(bx^c) \propto x^a Z_p(bx^c)$$
