You have the system
$$
\begin{align}
\dot{x}_1&=\sin(x_2)\\
\dot{x}_2&=-x_1-\sin(x_2)
\end{align}
$$
which has the linearization matrix
$$
A(x)=\begin{pmatrix}0&\cos{x_2}\\-1&-\cos{x_2}\end{pmatrix}
$$
which at the origin is
$$
A=A(0)=\begin{pmatrix}0&1\\-1&-1\end{pmatrix}
$$
Now taking $Q=I$ you can solve
$$
PA+A^TP=-Q
$$
for $P>0$ and get $P=\begin{pmatrix}3/2&1/2\\1/2&1\end{pmatrix}$ so that
$$
V(x)=x^TPx=3 x_1^2/2 + x_1 x_2 + x_2^2
$$
Now finding a region of attraction can be done by solving
$$
k=\min_x V(x) \text{ s.t. } \dot{V}(x)=0, x\neq0
$$
This is a constrained minimization problem:
$$
\begin{align}
k=&\min_x\quad 3 x_1^2/2 + x_1 x_2 + x_2^2 \\
&\text{ s.t. } \quad 2 x_1\sin(x_2) - 2 x_1 x_2 - x_2 \sin(x_2) - x_1^2=0\\
&\phantom{\text{ s.t. }}\quad x\neq0
\end{align}
$$
This is not easy to solve analytically but we can solve it numerically. We can solve $\dot{V}(x)=0$ for $x_2$ and get two different solutions which we can substitute in $V(x)$ to get:
$$
V_1(x_2)=\frac{\left(2\,x_{2}-3\,\sin\left(x_{2}\right)\right)\,\left(3\,x_{2}-2\,\sin\left(x_{2}\right)+2\,\sqrt{{x_{2}}^2-3\,x_{2}\,\sin\left(x_{2}\right)+{\sin\left(x_{2}\right)}^2}\right)}{2}
$$
and
$$
V_2(x_2)=-\frac{\left(2\,x_{2}-3\,\sin\left(x_{2}\right)\right)\,\left(2\,\sin\left(x_{2}\right)-3\,x_{2}+2\,\sqrt{{x_{2}}^2-3\,x_{2}\,\sin\left(x_{2}\right)+{\sin\left(x_{2}\right)}^2}\right)}{2}
$$
both assuming ${x_{2}}^2-3\,x_{2}\,\sin\left(x_{2}\right)+{\sin\left(x_{2}\right)}^2\geq 0.$ Now we can plot both functions:
In this image, only the parts of $V_1$ and $V_2$ are plotted where ${x_{2}}^2-3\,x_{2}\,\sin\left(x_{2}\right)+{\sin\left(x_{2}\right)}^2\geq 0.$ We can now see that we have two minima, which we can use to compute $x_1$. We end up with two solutions:
$$
\begin{align}
(x_1,x_2)&=(0.9885, -2.2013)\\
(x_1,x_2)&=(-0.9885, 2.2013)
\end{align}
$$
The first corresponds to the minimum of the blue graph, the second to the minimum of the red graph. In both cases we have
$$
k=V(0.9885,-2.2013)=V(-0.9885,2.2013)=4.1355
$$
which is the level you are looking for. We can confirm this by checking another plot:
- Teal: $\dot{V}(x)\leq 0$
- Dark blue: $\dot{V}(x)> 0$
- Yellow: $\dot{V}(x)\leq 0$ and $V(x)\leq k$
- Black cross: origin
- Blue dot: first solution $(x_1,x_2)=(0.9885, -2.2013)$
- Red dot: second solution $(x_1,x_2)=(-0.9885, 2.2013)$
Note however that the yellow area in general is only a subset of the region of attraction.
From the Jacobian that you obtained for the stable equilibrium points $(-1, 0)$, and $(1, 0)$, the linearized system becomes
$$ \left[\matrix{\dot{x}_{1} \cr \dot{x}_{2}}\right] = \left[\matrix{0 & 1 \cr -2 & -\delta}\right] \left[\matrix{x_{1} \cr x_{2}}\right] .$$
Consider the Lyapunov function in the following quadratic form:
$$ V(\mathbf{x}) = \frac{1}{2} \mathbf{x}^{T} \mathbf{P} \mathbf{x} = \frac{1}{2} \left[\matrix{x_{1} & x_{2}}\right] \left[\matrix{\frac{\delta^{2} + 2 \delta + 2}{\delta} & 1 \cr 1 & \frac{\delta + 1}{\delta}}\right] \left[\matrix{x_{1} \cr x_{2}}\right] .$$
The positive definite $\mathbf{P}$ is found by solving the Lyapunov equation $\mathbf{P} \mathbf{A} + \mathbf{A}^{T} \mathbf{P} = - \mathbf{Q}$. In this case, $\mathbf{Q}$ is chosen as
$$ \mathbf{Q} = \left[\matrix{2 & 0 \cr 0 & \delta}\right] .$$
The derivative of $V(\mathbf{x})$ along the trajectories of the linear system is given by
$$ \dot{V}(\mathbf{x}) = - 2 x_{1}^{2} - \delta x_{2}^{2} < 0 \quad \text{for} \quad \mathbf{x} \neq \mathbf{0} .$$
Remember that this result, however, is local.
Region of Attraction
To estimate the region of attraction, the previous Lyapunov function can be used to determine a domain $D$ about the stable equilibrium by solving the derivative of $V(\mathbf{x})$ along the trajectories of the nonlinear system (Duffing Oscillator) where $\dot{V}(\mathbf{x})$ is negative definite.
$$ \dot{V} = \left(\frac{\delta^{2} + 2 \delta + 2}{\delta}\right) x_{1} \dot{x}_{1} + \dot{x}_{1} x_{2} + x_{1} \dot{x}_{2} + \left(\frac{\delta + 1}{\delta}\right) x_{2} \dot{x}_{2} $$
$$ \dot{V} = \left(\frac{\delta^{2} + 2 \delta + 2}{\delta}\right) x_{1} x_{2} + x_{2}^{2} + x_{1} \left(x_{1} - x_{1}^{3} - \delta x_{2}\right) + \left(\frac{\delta + 1}{\delta}\right) x_{2} \left(x_{1} - x_{1}^{3} - \delta x_{2}\right) $$
$$ \dot{V} = - x_{1}^{4} - \left[\left(\frac{\delta + 1}{\delta}\right) x_{1} x_{2} - 1\right] x_{1}^{2} + \left(\frac{3 \delta + 3}{\delta}\right) x_{1} x_{2} - \delta x_{2}^{2} $$
For example, if $\delta = 1$, the region of attraction is found to be satisfied in $D_{(1, 0)} = \left\{\sqrt{2} < x_{1} < \sqrt{5}, x_{2} \in \mathbb{R}\right\}$, and $D_{(-1, 0)} = \left\{-\sqrt{5} < x_{1} < -\sqrt{2}, x_{2} \in \mathbb{R}\right\}$.
The streamplot for the Duffing Oscillator with $\delta = 1$ is shown below. Note that $\pm \sqrt{5} \approx \pm 2.236$.
Best Answer
Consider the same Lyapunov candidate, $$V(x) = \frac{1}{2} \left( x_1^2 + x_2^2 \right).$$
Differentiating along the vector field like you did,
$$ \begin{aligned} \dot{V}(x) &= -x_1^2+x_1\,x_2+x_1\,(x_1^2+x_2^2)\,\sin(t)-x_2\,x_1-x_2^2+x_2\,(x_1^2+x_2^2)\,\cos(t),\\ &= -x_1^2-x_2^2+x_1\,(x_1^2+x_2^2)\,\sin(t)+x_2\,(x_1^2+x_2^2)\,\cos(t)\end{aligned}$$
Now this is where I depart from your manipulation. Let me give you a hint to give you an idea of where I will take this.
Hint: $A \sin(t) + B \cos(t) = \sqrt{A^2 + B^2} \cos\left(t - \measuredangle A + j B\right)$ or, alternatively, use Cauchy-Scwharz as suggested by @Lutz (Sorry, I only saw your suggestion after I wrote this answer!).
Step 1:
Step 2: