Let us consider the heat IBVP:
\begin{align*}
u_t=\alpha^2u_{xx},\ 0<x<1,\ t>0\\
u_x(0,t)=0,\ t>0\\
u_x(1,t)+u(1,t)=0,\ t>0\\
u(x,0)=1,\ 0<x<1.
\end{align*}
First we solve the Sturm-Liouville problem:
\begin{align*}
X^{\prime\prime}=kX\ \mbox{with BCs}\\
\left\{
\begin{aligned}
&X'(0)=0\\
&X'(1)+X(1)=0.
\end{aligned}
\right.
\end{align*}
The two cases $k=0$ and $k=\lambda^2>0$ are left to be checked for students. Note that the case $k=0$ leads to a trivial solution and for $k=\lambda^2>0$ case there is no solution. Let $k=-\lambda^2$. Then
$$X(x)=A\cos\lambda x+B\sin \lambda x$$
and
$$X'(x)=-A\lambda\sin\lambda x+B\lambda\cos\lambda x.$$
We obtain $B=0$ from the condition $X'(0)=0$. So $X'(x)=-A\lambda\sin\lambda x$. The BC $X'(1)+X(1)=0$ then implies that the eigenvalues $\lambda$ satisfies
$$\lambda=\cot\lambda\ \ \ \ \ (1)$$

The roots of (1) can be found by the Newton’s method. In Maxima, the Newton’s method can be performed as follows.

First we define $\cot(x)-x$ as a function $f(x)$:

(%i1) f(x):=cot(x)-x;

Newton’s method applied to f(x) is performed by

(%i2) newton(t):=subst(t,x,x-f(x)/diff(f(x),x));

To use Newton’s method we need the initial approximation. We can make a guess from its graph. As you can see the first zero is near 0.5, so we choose $x_0=0.5$.

If you just run

(%i3) newton(0.5);

Maxima will return the output

(%o3)                              0.74865744241706

In order to calculate next approximation $x_1$, simply run

(%i4) newton(%);

(%o4)                              0.85239848545338

For the next approximation $x_2$, again run

(%i5) newton(%);

(%o5)                              0.86029880244504

By running the same command, the next approximation $x_3$ is obtained as $x_3=0.86033358835819$. This value is accurate to eight decimal places. If you think this is good enough, you can stop.

Similarly we obtain $\lambda_n$ for $n=2,3,4,\cdots$:
$$\begin{array}{|c||c|c|}
\hline
n & \lambda_n & \frac{(2n-1)\pi}{2}\\
\hline
1 & .8603358 & 1.570796327\\
\hline
2 & 3.425618 & 4.712388981\\
\hline
3 & 6.437298 & 7.853981635\\
\hline
4 & 9.529334 & 10.99557429\\
\hline
5 & 12.645287 & 14.13716694\\
\hline
\cdots & \cdots & \cdots\\
\hline
\end{array}$$
The table shows that as $n$ increases $\frac{(2n-1)\pi}{2}$ gets closer to $\lambda_n$. So for a practical purpose one may use the approximation
$$\lambda_n\approx \frac{(2n-1)\pi}{2},\ \lambda_n=1,2,3,\cdots$$
The temperature distribution $u(x,t)$ is then given by
$$u(x,t)=\sum_{n=1}^\infty A_ne^{-\alpha^2\lambda_n^2t}\cos\lambda_n x$$
The $L_n$’s and $A_n$’s are computed to be:
\begin{align*}
L_n&=\int_0^1\cos^2\lambda_nxdx\\
&=\frac{1}{2}\int_0^1[\cos(2\lambda_nx)+1]dx\\
&=\frac{1}{2}\left[\frac{\sin 2\lambda_n}{2\lambda_n}+1\right]\\
&=\frac{1}{2}\left[\frac{\sin\lambda_n\cos\lambda_n}{\lambda_n}+1\right]\ (\sin2\lambda_n=2\sin\lambda_n\cos\lambda_n)\\
&=\frac{1}{2}\left[\frac{\sin\lambda_n\cos\lambda_n+\lambda_n}{\lambda_n}\right]
\end{align*}
and
\begin{align*}
A_n&=\frac{1}{L_n}\int_0^1\phi(x)\cos\lambda_nxdx\\
&=\frac{2\lambda_n}{\sin\lambda_n\cos\lambda_n+\lambda_n}\int_0^1\cos\lambda_nxdx\\
&=\frac{2\sin\lambda_n}{\sin\lambda_n\cos\lambda_n+\lambda_n}.
\end{align*}

This calculation is done by assuming that $\cos\lambda_n x$, $n=1,2,\cdots$ are orthogonal functions. In fact they are orthogonal functions if $\lambda_n=\frac{(2n-1)\pi}{2}$, $n=1,2,\cdots$.

Finally we have
$$u(x,t)=\sum_{n=1}^\infty\frac{2\sin\lambda_n}{\sin\lambda_n\cos\lambda_n+\lambda_n}e^{-\alpha^2\lambda_n^2 t}\cos\lambda_n x.$$
For an approximation one may take
$$u(x,t)\approx\sum_{n=1}^\infty\frac{4(-1)^{n-1}}{(2n-1)\pi}e^{-\alpha^2\frac{(2n-1)^2\pi^2}{4}t}\cos\frac{(2n-1)\pi x}{2}\ \ \ \ \ (2)$$

The 3D plot of $u(x,t)$ in (2) is



References:

David Betounes, Partial Differential Equations for Computational Science with Maple and Vector Analysis, TELOS, Springer-Verlag