Solving differential equations of the form $\theta''(t)+\theta'(t)+\theta(t-\delta)=0$ TLDR:
Hi, is there an exam friendly* method to solve (without using numerical methods)  differential equations of the form
$$\theta''(t)+\theta'(t)+\theta(t-\delta)=0$$
I have been searching for 'delay differential equations' and only first order DDEs solutions seem to show up.
*possible with bare hands, a graphing calculator, (and a reasonable level of mathematical intelligence for a pretty smart high school student)

Basically, I am working on making a mock paper for my friends (still in high school) and ended up setting a question on 2nd order ODE that I can't solve. So I'm curious on whether there is a (hopefully not too advanced) way to solve such ODEs.
The full question and the solution I have worked on so far is below:

Consider an elliptical track at the standard position such that its major and minor axes lie on the $x$-axis and $y$-axis respectively.
A point $P\left(a\cos \theta,b\sin \theta\right)$ on the ellipse, where $a>b>0$, is moving with an angular acceleration $\alpha$ rad s-2, an angular velocity $\omega$ rad s-1, and an angular displacement $\theta$ rad measured in an anticlockwise manner from the starting point $P_0\left(a,0\right)$.
As $P$ is moving, a laser ray is continuously shot with the speed of light $c$ m s-1 from $P$ through the left focus of the elliptical track and reflects off a point $Q$ on the ellipse towards the right focus. The right focus then directs the laser ray back towards the point where it was shot from. The time taken in the process is denoted by $\delta$ s. The angular displacement of the current position of $P$ with respect to the point where the laser eventually hits is $\phi$ rad.
Given that
$$\alpha+\omega+\theta=\phi$$
for time $t\ge\delta$, and 
$$\theta=\frac{2\sqrt3}3e^{-\frac t2}\sin\frac{\sqrt3}2t$$
for $0\le t <\delta$, find the parametric equations describing the motion of $P$ with respect to time $t$. 

By the geometric definition of an ellipse,
$PF_1+PF_2=2a$.
Hence, the distance traveled by the laser is 
\begin{align*}
PF_1+F_1Q+QF_2+F_2P
&=\left(PF_1+PF_2\right)+\left(QF_1+QF_2\right)\\
&=4a
\end{align*}
The time $\delta$ s taken for the laser to travel back is $$\delta=\frac{4a}{c}$$
Hence the angular distance is such that
$$\phi=\theta\left(t\right)-\theta\left(t-\delta\right)=\theta\left(t\right)-\theta\left(t-\frac{4a}{c}\right)$$
Rewriting the given differential equation,
\begin{align*}
\alpha+\omega+\theta&=\phi\\
\alpha\left(t\right)+\omega\left(t\right)
+\theta\left(t\right)&=\theta\left(t\right)-\theta\left(t-\delta\right)\\
\alpha\left(t\right)+\omega\left(t\right)
+\theta\left(t-\delta\right)&=0\\
\alpha\left(t-\delta\right)+\omega\left(t-\delta\right)
+\theta\left(t-2\delta\right)&=0\\
\end{align*}

 A: Using the Laplace transform we have
$$
s^2\Theta(s)+s\Theta(s)+e^{-\delta s}\Theta(s) = \dot\theta_0 + \theta_0+s\theta_0
$$
or
$$
\Theta(s) = \frac{\dot\theta_0 + \theta_0+s\theta_0}{s^2+s+e^{-\delta s}}
$$
and now taking a Padé expansion for $e^{-\delta s}$
NOTE
$$
e^{-\tau s} = \frac{-\frac{1}{120} s^3 \tau ^3+\frac{s^2 \tau ^2}{10}-\frac{s \tau }{2}+1}{\frac{s^3 \tau ^3}{120}+\frac{s^2 \tau ^2}{10}+\frac{s \tau }{2}+1}
$$
is a $3$th order Padé expansion.
Attached a plot for $\delta = \{0, 0.1,0.2,0.3,0.4,0.5,0.6\}$ The red one is for $\delta = 0$

NOTE
In the plot $\theta_0 = 1, \dot\theta_0 = 0$
A: The equations you set up can be summed up as.
\begin{cases} \theta(0) = 0,\ \dot{\theta}(0) = 1 \\ \ddot{\theta}(t) + \dot{\theta}(t) + \theta(t) = 0, & 0 \le t < \delta \\ \ddot{\theta}(t) + \dot{\theta}(t) + \theta(t-\delta) = 0, & t \ge \delta \end{cases}
A common tactic for delayed ODE is to solve it "piece-wise". You already have the solution in $[0,\delta)$, which I'll denote as
$$ \theta_0(t) = \frac{2}{\sqrt{3}}e^{-t/2}\sin\left(\frac{\sqrt{3}t}{2}\right) $$
You can then find a solution in $[\delta,2\delta)$, which I'll denote as $\theta_1(t)$. Since $\theta(t-\delta)_{t\in[\delta,2\delta)} = \theta_0(t-\delta)$ is a known function, the problem is reduced down to a regular, inhomogeneous ODE
$$ \ddot\theta_1(t) + \dot\theta_1(t) =-\theta_0(t-\delta), \qquad \theta_1(\delta) = \theta_0(\delta), \ \dot\theta_1(\delta) = \dot\theta_0(\delta) $$
This is in fact a first-order ODE in $\dot\theta_1$, so you have
$$ \frac{d}{dt}(e^t\dot\theta_1) = -e^t\theta_0 $$
$$ e^t\dot\theta_1 = -\int_\delta^t e^u\theta_0(u-\delta)\ du+C $$
Given $\dot\theta_1(\delta) = \theta_0(\delta) $, we get $C = \dot\theta_0(\delta)e^\delta$, therefore
$$ \dot\theta_1(t) = \dot\theta_0(\delta)e^{-(t-\delta)} - e^{-t}\int_{\delta}^t e^u\theta_0(u-\delta)\ du $$
Integrating once more and using $\theta_1(\delta) = \theta_0(\delta)$
\begin{align} \theta_1(t) &= \theta_0(\delta) + \int_{\delta}^t \dot\theta_1(s)\ ds \\
&= \theta_0(\delta) + \dot\theta_0(\delta)\big[1-e^{-(t-\delta)}\big] - \int_{\delta}^te^{-s}\int_{\delta}^s e^u\theta_0(u-\delta)\ du\ ds 
\end{align}
In fact, all solutions in successive intervals have this form. If we denote $\theta_n(t)$ as the solution in the interval $[n\delta, (n+1)\delta)$, then the following recurring relation holds
$$ \theta_n(t) = \theta_{n-1}(n\delta) + \dot\theta_{n-1}(n\delta)\big[1-e^{-(t-n\delta)}\big] - \int_{n\delta}^t e^{-s}\int_{n\delta}^s e^u\theta_{n-1}(u-\delta)\ du\ ds $$
The integrals are doable by hand, but they will get ugly.
