Issues carrying out Taylor Expansion for small oscillations I have an equation of motion for a pendulum: $\ddot{\theta}+\frac{g}{l}\sin(\theta)-\dot{\phi}^2\sin(\theta)\cos(\theta)=0$
I want to taylor expand this equation for small oscillations, i.e.: $\theta = \theta_0 + \delta\theta$, and I should end up with the following expression: $\delta\ddot{\theta}+\dot{\phi}_0^2(1+3\cos^2(\theta_0))\delta\theta \approx 0$ where $\dot{\phi}_0 = \sqrt{g/(l\cos(\theta_0))}$ is constant
I've tried to a number of times, but I can't seem to get this to work at all, for simpler taylor expansions I don't have any issue but this one is really stumping me.
I appreciate any help, and thanks in advance.
ETA: My apologies, here's some more information:
I'm dealing with a spherical pendulum.
$\theta$ is the angle the (massless) rod of length $l$ holding the mass makes with the $z$-axis.
$\phi$ is the angle the mass makes in the $x$-$y$ plane with regards to the positive $x$-axis.
$g$ is gravity
 A: I will first describe how I understand the question: You have a differential equation, namely $\ddot{\theta}+\left(\frac{g}{l}-\dot{\phi}^2\cos(\theta)\right)\sin(\theta)=0$. You mention that you want to work out what happens for small oscillations, so you probably want to expand around the equilibrium position. "Equilibrium" means that the second time derivative of the function we are dealing with is zero, so $\ddot{\theta}=0$. Plugging in our differential equation we get $\frac{g}{l}-\dot{\phi_0}^2\cos(\theta_0)=0$ where $\theta_0$ and $\phi_0$ indicate the values of $\theta$ and $\phi$ in the equilibrium position (around which you wish to expand). This finally yields $\dot{\phi_0}^2=\left(\frac{g}{l}\cos(\theta_0)\right)^{\frac{1}{2}}$. This is exactly what you provided so I assume that I got everything right until now. From now on I try to do the calculation while only skipping minor steps (like multiplying stuff out) and only using the Taylor series os $\sin$ and $\cos$ and one additional information that was not provided in your question but that has to be used (I will say what it is later on).
We plug in $\theta=\theta_0+\delta\theta$. Because one now has a sum inside the trigonometric functions it is useful to consider the trigonometric identities $\sin(\alpha+\beta)=\sin(\alpha)\cos(\beta)+\cos(\alpha)\sin(\beta)$ and $\cos(\alpha+\beta)=\cos(\alpha)\cos(\beta)-\sin(\alpha)\sin(\beta)$ with $\alpha=\theta_0$ and $\beta=\delta\theta$.
Our differential equation becomes $$\delta\ddot{\theta}+\left(\frac{g}{l}-\dot{\phi}^2\cos(\theta_0)\cos(\delta\theta)+\dot{\phi}^2\sin(\theta_0)\sin(\delta\theta)\right)(\sin(\theta_0)\cos(\delta\theta)+\cos(\theta_0)\sin(\delta\theta))=0$$ where I used $\ddot{\theta}_0=0$ because it is a constant.
Until now everything was exact and no approximation has taken place. Now we want to approximate the functions containing $\delta\theta$ by some Taylor polynomial. The first relevant terms of the Talor series for $\sin$ and $\cos$ around zero should be known. Why do we now suddenly expand around zero? Above we rewrote the differential equations in a way to separate $\theta_0$ from $\delta\theta$ so now we are able to look at the functions containing only $\delta\theta$.
The known expansions are: $\cos(\delta\theta)=1-\frac{(\delta\theta)^2}{2}+...$ and $\sin(\delta\theta)=\delta\theta+...$ We will only take the leading term so actually we can ignore the quadratic term in the first one: $\cos(\delta\theta)=1+...$ where I used "$...$" to denote higher order terms which we do not care about in our approximation.
But we have to be careful: $\dot{\phi}^2$ also has to be expanded. I have intentionally ignored that point until now. Without that additional fact one gets to a "wrong" conclusion, like e.g. user @Eli in another answer (which is not their fault because they just expanded the parts of the differential equation given in your question while one needs additional information). The additional information we have is that the angular momentum $L=\sin^2(\theta)\dot{\phi}$ (I am ignoring come factors here which are not important) is constant along a trajectory. We can assume that at some point in the trajectory the system passes through the equilibrium, i.e. we have a point in time where $\theta=\theta_0$ and $\dot{\phi}=\dot{\phi}_0$. This means the conservation of angular momentum takes the form $\sin^2(\theta)\dot{\phi}=\sin^2(\theta_0)\dot{\phi}_0$. We know already that $\sin(\theta)=\sin(\theta_0)+\cos(\theta_0)\delta\theta+...$, so it follows from the binomial theorem that $\sin^2(\theta)=\sin^2(\theta_0)+2\sin(\theta_0)\cos(\theta_0)\delta\theta+...$ (where we ignored the term quadratic in $\delta\theta$). Also we have $\dot{\phi}=\dot{\phi}_0+\delta \dot{\phi}$ where $\delta\dot{\phi}$ is actually the thing we need because we want to expand $\dot{\phi}$ around the equilibrium position. Multiplying everything out, again leaving out the terms quadratic in $\delta\theta$, we obtain $\delta\dot{\phi}=-2\frac{\cos(\theta_0)}{\sin(\theta_0)}\dot{\phi}_0\delta\theta$.
To include that in our differential equation we again use the binomial theorem to get $\dot{\phi}^2=\dot{\phi}_0^2+2\dot{\phi}_0\delta\dot{\phi}+...=\dot{\phi}_0^2-4\frac{\cos(\theta_0)}{\sin(\theta_0)}\dot{\phi}_0^2\delta\theta+...$
Plugging that in yields $$\delta\ddot{\theta}+\left(\frac{g}{l}-\left(\dot{\phi}_0^2-4\frac{\cos(\theta_0)}{\sin(\theta_0)}\dot{\phi}_0^2\delta\theta\right)(\cos(\theta_0)-\sin(\theta_0)\delta\theta)\right)(\sin(\theta_0)+\cos(\theta_0)\delta\theta)\approx 0$$ When multiplying that out we remember that we will only use terms that are exactly linear in $\delta\theta$. The constant terms will cancel because we did expand around the minimum and everything quadratic or higher will be ignored. This gives finally $$\delta\ddot{\theta}+\left(\frac{g}{l}\cos(\theta_0)+\dot{\phi}_0^2 (3\cos(\theta_0)^2+\sin^2(\theta_0))\right)\delta\theta \approx 0$$. This is almost what you want. Now only insert the equilibirum condition $\frac{g}{l}\cos(\theta_0)=\dot{\phi}_0^2 \cos(\theta_0)^2=\dot{\phi}_0^2 (1-\sin(\theta_0)^2)$ (the last step is Pythagoras). With that we have the desired $\delta\ddot{\theta}+\dot{\phi}_0^2(1+3\cos^2(\theta_0))\delta\theta\approx 0$. I hope that was helpful.
A: The nonlinear function you want to linearize is:
$$f(\theta,\ddot{\theta},\dot{\phi}) = \ddot{\theta} + \frac{g}{l}\sin(\theta) - \dot{\phi}^2 \sin(\theta)\cos(\theta) = 0 \tag 1$$
Linearization is done by expanding the function in Taylor series but only up to the first term:
$$f_L(\theta,\ddot{\theta},\dot{\phi}) = f(\theta_0,\ddot{\theta}_0,\dot{\phi}_0) + \left.\frac{\partial}{\partial\theta}f(\theta,\ddot{\theta},\dot{\phi})\right|_{P_0} \Delta\theta + \left.\frac{\partial}{\partial\ddot{\theta}}f(\theta,\ddot{\theta},\dot{\phi})\right|_{P_0} \Delta\ddot{\theta} + \left.\frac{\partial}{\partial\dot{\phi}}f(\theta,\ddot{\theta},\dot{\phi})\right|_{P_0} \Delta\dot{\phi}$$
where $P_0 = (\theta_0,\ddot{\theta}_0,\dot{\phi}_0)$ is the equilibrium point around which linearization is performed, and $f(\theta_0,\ddot{\theta}_0,\dot{\phi}_0) = 0$ is obvious from Eq. (1). By solving the above equation we obtain:
$$f_L(\theta,\ddot{\theta},\dot{\phi}) = \Delta\ddot{\theta} + \left(\frac{g}{l}\cos(\theta_0) - \dot{\phi}^2_0 \left(\cos^2(\theta_0) - \sin^2(\theta_0)\right) \right) \Delta\theta + \left(-2\dot{\phi}_0\sin(\theta_0)\cos(\theta_0) \right) \Delta \dot{\phi}$$
where the equilibrium point is defined as:
$$\ddot{\theta}_0 = 0, \qquad \dot{\phi}^2_0 = \frac{g}{l}\frac{1}{\cos(\theta_0)}$$
Now plug $\dot{\phi}_0^2$ in the linearized equation and you get:
$$f_L(\theta,\ddot{\theta},\dot{\phi}) = \Delta\ddot{\theta} + \dot{\phi}_0^2 \left(\frac{g}{l}\cos(\theta_0)\frac{l}{g}\cos(\theta_0) - \left(\cos^2(\theta_0) - \sin^2(\theta_0)\right)\right) \Delta\theta + \left(-2\dot{\phi}_0\sin(\theta_0)\cos(\theta_0) \right) \Delta\dot{\phi}$$
Finally, the linearized differential equation is
$$f_L(\theta,\ddot{\theta},\dot{\phi}) = \Delta\ddot{\theta} + \left(\dot{\phi}_0^2 \sin^2(\theta_0) \right)\Delta\theta + \left(-2\dot{\phi}_0\sin(\theta_0)\cos(\theta_0) \right) \Delta\dot{\phi} = 0 \tag 2$$
The above equation does not include any additional information about the $\dot{\phi}$.

As pointed out by @unsure the angular momentum is constant. This is second equation of the model that can be linearized:
$$h(\theta,\dot{\phi}) = \sin^2(\theta) \dot{\phi} + \text{const.} = 0 \tag 3$$
Following the same linearization approach:
$$h_L(\theta,\dot{\phi}) = h(\theta_0,\dot{\phi}_0) + \left.\frac{\partial}{\partial\theta}h(\theta,\dot{\phi})\right|_{P_0} \Delta\theta + \left.\frac{\partial}{\partial\dot{\phi}}h(\theta,\dot{\phi})\right|_{P_0} \Delta\dot{\phi}$$
where $h(\theta_0,\dot{\phi}_0) = 0$ is obvious from Eq. (3), we obtain:
$$h_L(\theta,\dot{\phi}) = \left(2\sin(\theta_0)\cos(\theta_0)\dot{\phi}_0\right) \Delta\theta + \left(\sin^2(\theta_0)\right) \Delta\dot{\phi} = 0$$
from which the following identity is obtained:
$$\Delta\dot{\phi} = \left(-2\frac{\cos(\theta_0)}{\sin(\theta_0)}\dot{\phi}_0\right) \Delta\theta \tag 4$$

By plugging Eq. (4) into Eq. (2) and after doing some simple arithmetic we obtain:
$$\Delta\ddot{\theta} + \left(\dot{\phi}^2_0 \left(1 + 3\cos^2(\theta_0) \right)\right) \Delta\theta = 0 \quad \surd$$
This equals the form given in the original post.
There are different ways to linearize nonlinear models, but in my experience expanding the function in Taylor series up to the first term is the simplest and most straightforward.
A: this is the function that you want to linearized
$$f={\frac {g\sin \left( \theta  \right) }{l}}-\frac 12\,{\dot\phi }^{2}\,\sin(\theta)\,\cos(\theta)
\tag 1$$
from the conservation of the energy $~E=E_0~$ you obtain $~\dot{\phi}^2~$
where
\begin{align*}
&E= \frac m2\left(  \left( \sin \left( \theta \right)  \right) ^{2}{\dot\phi }^{2}+{
\dot\theta }^{2} \right) {l}^{2}-m\,g\,l\,\cos(\theta)
\\
&E_0=E(\theta=\theta_0~,\dot{\phi}=\dot{\phi}_0)\quad\Rightarrow\\
&\dot{\phi}^2={\frac { \left( \sin \left( \theta_{{0}} \right)  \right) ^{2}{\dot\phi _
{{0}}}^{2}}{ \left( \sin \left( \theta \right)  \right) ^{2}}}-2\,{
\frac {g \left( -\cos \left( \theta \right) +\cos \left( \theta_{{0}}
 \right)  \right) }{l \left( \sin \left( \theta \right)  \right) ^{2}}
}
\end{align*}
substitute into equation  (1)
with $~\theta=\theta_0+\delta\theta~$ and linearized the result    you obtain
\begin{align*}
 &f_L={\frac {\delta \theta \, \left( l{\dot\phi _{{0}}}^{2}+3\,g\cos \left(
\theta_{{0}} \right)  \right) }{l}}-\sin \left( \theta_{{0}} \right)
\cos \left( \theta_{{0}} \right) {\dot\phi _{{0}}}^{2}+{\frac {g\sin
 \left( \theta_{{0}} \right) }{l}}
\end{align*}
With
$~f_L=0\quad,\delta\theta=0\quad\Rightarrow~$ $~\dot{\phi}_0^2=\frac{g}{\cos(\theta_0)\,l}~$ hence
\begin{align*}
  &f_L\left(\dot{\phi}_0^2=\frac{g}{\cos(\theta_0)\,l}\right)={\frac {\delta \theta \,g \left( 3\, \left( \cos \left( \theta_{{0}}
 \right)  \right) ^{2}+1 \right) }{\cos \left( \theta_{{0}} \right) l}
}=\dot{\phi}_0^2\,\delta\theta\,\left(3\,\cos(\theta_0)^2+1\right)\\
&\delta\ddot\theta+\dot{\phi}_0^2\,\delta\theta\,\left(3\,\cos(\theta_0)^2+1\right)=0\quad\surd
 \end{align*}

The energy and equations of motoin
\begin{align*}
 &\text{the postion vector to the mass}\\
 &\mathbf{R}=l\,\left[ \begin {array}{c} \cos \left( \phi \right) \sin \left( \theta
 \right) \\  \sin \left( \phi \right) \sin \left(
\theta \right) \\  \cos \left( \theta \right)
\end {array} \right]
 \end{align*}
from here you obtain the velocity
\begin{align*}
   &v_i =\frac{\partial R_{i}}{\partial q_j}\,\dot{q}^j\quad,
   \text{where}
   \quad\mathbf{q}=\begin{bmatrix}
                 \phi \\
                 \theta \\
               \end{bmatrix}
\quad\mathbf{\dot{q}}=\begin{bmatrix}
                 \dot\phi \\
                 \dot\theta \\
               \end{bmatrix}   \quad
\text{hence, the kinetic energy}\\
&T=\frac{m}{2}\,v_i\,v^i=
\frac{m\,l^2}{2}\left(\,\sin^2(\theta)\,\dot{\phi}^2
+\dot{\theta}^2\right)\quad \text{the potential energy}\\
&U=-m\,g\,l\,\cos(\theta)\quad \text{and the total energy}\\\\
&E=T+U=\frac{m\,l^2}{2}\left(\,\sin^2(\theta)\,\dot{\phi}^2\right)
+\dot{\theta}^2\,-m\,g\,l\,\cos(\theta) \\                          
 \end{align*}
the equations of motion with EL
\begin{align*}
  &\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{q}_i}\right)+\frac{\partial U}{\partial q_i}=0 \quad\Rightarrow\\\\
  &\ddot\theta -\sin \left( \theta \right) \cos \left( \theta \right)
{\dot\phi }^{2}+{\frac {\sin \left( \theta \right) g}{l}}
=0\\
&\ddot\phi +2\,{\frac {\cos \left( \theta \right) \dot\theta \,\dot\phi }{
\sin \left( \theta \right) }}=0\\\\
&\text{$~\phi~$ is cyclic coordinate}\quad \Rightarrow\\
&\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{\phi}}\right)=0\quad\Rightarrow\quad \sin^2(\theta)\,\dot{\phi}=\text{constant}
 \end{align*}
