There is a first error in your first equation
$$
\ddot{\pmb r}(t)=\left(-l\sin\theta\;\dot\theta^2+l\cos\theta\;\ddot\theta,\;l\cos\theta\;\dot\theta^2+l\sin\theta\;\ddot\theta\right)
$$
that is $$
m\ddot{\pmb r}(t)=-mg\hat y+\pmb T
$$
and we have
\begin{align}
ml(\ddot{\theta}\cos\theta-\dot{\theta}^2\sin\theta)&=-T\sin\theta. & \tag {eq. 1} \\
ml(\ddot{\theta}\sin\theta+\dot{\theta}^2\cos\theta )&=T \cos\theta - mg & \tag {eq. 2}
\end{align}
Observe that $T\neq mg\cos\theta$ (this is your second error).
Taking the difference $\sin\theta\times (\text{eq. }1)-\cos\theta \times (\text{eq. }2)$ we find
$$ ml\dot\theta^2=T-mg\cos\theta \tag 3$$
and taking the sum $\cos\theta\times (\text{eq. }1)+\sin\theta \times (\text{eq. }2)$ we find
$$ ml\ddot\theta=-mg\sin\theta \tag 4$$
that is
\begin{align}
ml\dot\theta^2+mg\cos\theta &=T\tag 5\\
\ddot\theta+\frac{g}{l}\sin\theta &=0 \tag 6
\end{align}
I think you can solve these equations.
You already laid out that the acceleration along the path has to be tangent to the path, $m\vec a=(\vec F\vcenter{\huge⋅}\hat u)\,\hat u$, as all force components perpendicular to the path are counter-acted by the path itself, be it as physical construction or force field.
As $\vec a$ is also the second time derivative of $\vec q(θ(t))$, one gets
$\vec a=\vec q''(θ)\dot θ^2+\vec q'(θ)\ddot θ$, where here $(\vec q''(θ)\vcenter{\huge⋅}\vec q'(θ))=0$. For the projection we thus get
$$
(\vec a \vcenter{\huge⋅}\vec q'(θ))=\|\vec q'(θ)\|^2\ddot θ
\\
=(\vec F/m \vcenter{\huge⋅}\vec q'(θ))=-g·R\cos(θ)
$$
Using $\|\vec q'(θ)\|^2=R^2+k^2$, this results in
$$
\ddotθ= -\frac{g \cdot R \cdot \cos(\theta)}{R^2 + k^2}
$$ along the tangent direction.
So what this reduces to is a pendulum equation with the increased length $\frac{R^2+k^2}R$, the spiral being oriented horizontally, so the system resembles a loop-de-loop. Note that $θ=0$ is the horizontal $x$ direction.
To get a function table in $(t,x,y,z)$ for the motion, you need to integrate this second order equation and then insert the function table for $θ(t)$ into the formula for the cylinder coordinates.
One could, as you did, add the derivatives of the cylinder coordinates to the differential system. But using the low-order Euler method would introduce additional errors, moving away from the spiral, over those from the angle integration.
As it is a Hamiltonian system, you can easily change the Euler method to semi-implicit or symplectic Euler, or with another small twist, into the leapfrog Verlet method.
Best Answer
Mathematically, not really.
I noticed that your formula for pendulum is not correct. You might mean these: $$\frac{\dot{\theta}^2}{2} = \frac{g}{l} \cos(\theta)$$ and $$\ddot{\theta} = -\frac{g}{l} \sin(\theta)$$
The trivial counterexample is $\theta \equiv \pi /2$, where $\dot{\theta}^2/2=\frac{g}{l}\cos\theta = 0$ but $\ddot{\theta}=0 \ne -\frac{g}{l}\sin\theta$.
In that case, you may assume that $\theta$ is not a constant.
By taking derivative, you get $\ddot{\theta}\dot{\theta}=\frac{g}{l}\sin\theta \dot{\theta}$ and you can cancel $\dot{\theta}$ because $\theta$ is not a constant, then you get $\ddot{\theta} = -\frac{g}{l}\sin\theta$. However, if you have $\ddot{\theta} = -\frac{g}{l}\sin\theta$, then you can't get $\frac{\dot{\theta}^2}{2} = \frac{g}{l} \cos(\theta)$, because actually $\frac{\dot{\theta}^2}{2} = \frac{g}{l} \cos(\theta) + C$.