Tension in a simple pendulum

You are doing nothing wrong - there is no simple exact solution to this problem.

However it's easily solvable if you use the small angle approximation,

$$ sin(\theta) \approx \theta $$ $$ cos(\theta) \approx 1 - \frac{1}{2}\theta^2$$

$\ddot{\theta} = -\frac{g}{l}\theta$ has an easy solution,

$$\theta(t) = A\sin(\omega t)$$

where $\omega = \sqrt{\frac{g}{l}}$

It helps me to remember sin and cos are just infinite series!

$$ sin(x) = x - \frac{1}{3!}x^3 + \frac{1}{5!}x^5 - \frac{1}{7!}x^7 + ...$$

$$ cos(x) = 1 - \frac{1}{2!}x^2 + \frac{1}{4!}x^4 - \frac{1}{6!}x^6 +...$$

So when x is small these approximations make sense.

There are two reasons that your equations came out differently that you expected. One is a simple geometry error: you should have $T_x = - T \sin \theta$ and $T_y = - T \cos \theta$. The other is that you used Cartesian coordinates to write down Newton's Second Law rather than polar (which is the usual choice.) But it's not too hard to take your equations (corrected appropriately) and show that they imply the more customary form.
\begin{align} m(l\ddot{\theta}\cos\theta-l\dot{\theta}^2\sin\theta)&=-T\sin\theta. & (\text{Eq}~1) \\ m(-l\ddot{\theta}\sin\theta-l\dot{\theta}^2\cos\theta )&=-T \cos\theta + mg & (\text{Eq}~2) \end{align} Taking the sum of $(\sin \theta) (\text{Eq}~1) + (\cos \theta) (\text{Eq}~2)$, we obtain $$ m(l\ddot{\theta}\cos\theta\sin\theta -l\dot{\theta}^2\sin^2\theta)+ m(-l\ddot{\theta}\sin\theta\cos\theta-l\dot{\theta}^2\cos^2\theta )=-T\sin^2\theta-T \cos^2\theta + mg \cos \theta \\ \Rightarrow - m l \dot{\theta}^2 = -T + mg \cos \theta. $$ Similarly, taking the sum of $(\cos \theta) (\text{Eq}~1) - (\sin \theta) (\text{Eq}~2)$, we get $$ m(l\ddot{\theta}\cos^2\theta -l\dot{\theta}^2\sin\theta \cos \theta)- m(-l\ddot{\theta}\sin^2\theta-l\dot{\theta}^2\cos\theta \sin\theta )=-T\sin\theta \cos \theta + T \cos\theta \sin \theta - mg \sin \theta \\ \Rightarrow m l \ddot{\theta} = - mg \sin \theta. $$

This is an illustration of a general principle, by the way: the equations of motion expressed in one set of coordinates will always be equivalent to the equations of motion expressed in another set of coordinates. In particular, in your case you were expressing the equation $\vec{F} = m \vec{a}$ in terms of its Cartesian components ($F_x = m a_x $ and $F_y = m a_y$), while the usual form is in polar coordinates. But the components $\vec{F}$ in polar coordinates are related to the components in Cartesian components by $$ \begin{bmatrix} F_r \\ F_\theta \end{bmatrix} = \begin{bmatrix} \sin \theta & \cos \theta \\ - \cos \theta & \sin \theta \end{bmatrix} \begin{bmatrix} F_x \\ F_y \end{bmatrix} $$ and similarly for $\vec{a}$. Thus, by taking the appropriate linear combinations of your equations of motion, we can express them in polar coordinates instead.

For an undamped pendulum, we can actually have an exact, but not closed-form solution but this is better than approximating the model, but the simple pendulum often approximates small angles. I could also say the simple pendulum is simple because it is undamped, or because it is a point mass in plane rotation. Before I go into that, here are the derivations of the solutions you were looking for. Also, I hope you can tell me what software you use to diagram, as it looks much better than what I usually do.

Using the $\theta$ coordinate you specified, we can formulate the equation of motion from the torque equation about point A:

$$ I^A\ddot{\theta} = \sum \vec{T^A} \\ I^o\ddot{\theta} = -mgl~sin\theta $$

But for a point mass, we can say $I^o=\int_Mr^2~dm=ml^2$ and this can take us to one of the answers you were looking for.

$$ ml^2\ddot{\theta} = -mgl~sin\theta \\ \ddot{\theta} = -\frac{g}{l}sin\theta $$

If we wanted to get the other answer, then we can look at the acceleration of the system and use newton's 2nd Law for a point mass. $\vec{F}=m\vec{a}$. Where the ball must accelerate to navigate the turn at an acceleration that can be found to be for any turn $a_c=\omega^2R$. In this situation, $R=l$ and $\omega=\dot{\theta}$. Writing the equation in an axis that align with the tension,

$$ ma = ma_c = T - mg~cos\theta \\ T = ma_c+mg~cos\theta = m\dot{\theta}^2l + mg~cos\theta $$

Okay, to get to the solution for the equation of motion we don't need tension since this equation has enough to describe the system, further, since tension acts perpendicular to displacement, it contributes no work and is not required to be used to calculate the state of the system.

$$ \ddot{\theta} + \frac{g}{l}sin\theta = 0 $$

By inspiration, or staring at it long enough you might decide to do this, innocently, not so innocently, who knows?

$$ \dot{\theta}\ddot{\theta} + \dot{\theta}\frac{g}{l}sin\theta = 0 $$

It doesn't seem easier, but from here you could say that the above is the same as the first thing below by the chain rule. From here I'll just work it to the end, I'll call the constant from integration $C$.

$$ \frac{d}{dt} \bigg[\dot{\theta}^2/2 - \frac{g}{l}cos\theta\bigg] = 0 \\ \dot{\theta}^2/2 - \frac{g}{l}cos\theta = C \\ \frac{d\theta}{dt} = \sqrt{ 2C+2\frac{\strut g}{l}cos\theta} \\ \int d\theta\bigg/\sqrt{ 2C+2\frac{\strut g}{l}cos\theta} = \int{dt} $$

So it might be nasty to solve, but this isn't based on an approximation, so the model is better. In practice, you'll find solving either of the ODE's the first order or the second order (the original) to be faster than computing this for $\Delta\theta$'s but wolfram gave a solution to it with an elliptic integral solution, so maybe it's alright, but it would still be $t = f(\theta)$ where we would be unable to construct an $f^{-1}$ such that $\theta = f^{-1}(x)$.

The last thing to do is to discuss if that integration constant has any meaning. If I'm to be honest, from this point I will be guessing/trying to remember, but we found that there was a value that was constant. The first line of the last series of equations can be written as follows if we use the second line:

$$ \frac{d}{dt}~C = 0 $$

So we know that it is a conserved quantity. Since it has to do with the square of speed I would expect that $C$ is proportional to energy rather than momentum, since momentum is linear with speed. This also has a good reason to be energy, if the energy is too small, the integral is not defined for all values of $\theta$ which means you won't have enough energy to get to certain heights.There exists a critical value of C such that we can get exactly to $\theta=\pi$, which is $C=g/l$ and this would correspond to the energy required to be at or attain unstable equilibrium - where the pendulum is balanced up at the top, with no way to fall down and the integrand becomes undefined for that value of $\theta$.