In principle you could just integrate the coupled Newton-Lorentz-Maxwell equations, but this is hard. :) If radiation reaction isn't important then this is just the regular two body problem with the Coulomb force rather than the gravitational force. In actual practice I would use a perturbation theory based on the two body problem with the radiation reaction force as a small perturbation.
You could develop a systematic perturbation theory if you want. The key idea is well described by Wikipedia:
Perturbation theory comprises mathematical methods that are used to find an approximate solution to a problem which cannot be solved exactly, by starting from the exact solution of a related problem. Perturbation theory is applicable if the problem at hand can be formulated by adding a "small" term to the mathematical description of the exactly solvable problem.
Here the related problem with an exact solution is the motion of a free particle and the small parameter is the small acceleration experienced by the particle. What you would do is expand every quantity in a series of ever small corrections. We can use a parameter $\lambda$ to count the order of smallness (at the end of the day we'll set $\lambda=1$):
$$\begin{array}{rl}
r &= r_0 + \lambda r_1 + \lambda^2 r_2 + \cdots \\
v &= v_0 + \lambda v_1 + \lambda^2 v_2 + \cdots \\
a &= \lambda a_1 + \lambda^2 a_2 + \cdots \\
\end{array}$$
The key thing to notice is that the acceleration starts out with an order of smallness. You'll substitute these into the exact equations of motion, expand everything out, and match term of the same order. What you'll get is a series of equations: the first one describing the motion of a free particle, the second one describing the first order correction due to the acceleration, the third one describing the next order correction due to the fact that the corrected trajectory changes the acceleration slightly (since $r_1$ feeds into $a_2$), and so on... You have presumably already solved the first two equations of this hierarchy.
This particular example is a little messy. In case you haven't seen any of this before I'll illustrate the basic mechanics with a simple example where we know the exact solution. Consider the quadratic equation
$$ x^2 - 2\epsilon x - 1 =0, $$
where $\epsilon$ is a small real number representing the perturbation. We know the exact solutions are
$$ x = \epsilon \pm \sqrt{1+\epsilon^2}. $$
But suppose we didn't know that. We could still make progress by writing
$$ x = x_0 + \epsilon x_1 + \epsilon^2 x_2 + \cdots. $$
Substituting that into the quadratic equation and gathering terms by powers of $\epsilon$ gives
$$ (x_0^2-1) + \epsilon (2 x_0 x_1 - 2 x_0) + \epsilon^2 (2 x_0 x_2 + x_1^2 - 2 x_1) + \cdots = 0, $$
where we only keep terms up to $\epsilon^2$ because we only expanded $x$ that far. Now every order in $\epsilon$ has to vanish seperately. (Why? Because if you had a cancellation between terms of different order that would mean your coefficients are vastly different sizes, meaning your perturbation is not really having a small effect on the solution, i.e. the perturbation theory is fundamentally broken. This is the case for a structurally unstable system: a small perturbation changes the behaviour completely. You could try identifying another small parameter where the scaling is better behaved.) So this gives you a series of equations:
$$\begin{array}{rl}
x_0^2 - 1 &= 0, \\
2 x_0 x_1 - 2 x_0 &= 0, \\
2 x_0 x_2 + x_1^2 - 2 x_1 &= 0, \\
\vdots
\end{array}$$
The first equation is the unperturbed problem with the exact solutions $x_0 = \pm 1$. Note that the rest of the equations determine the next order of the solution in terms of the lower order quantities. Now there are two choices of $x_0$ to base the perturbation series on. In this case they both give valid solutions. I'll pick $x_0 = +1$ for now and let you develop the other solution. Now we use this $x_0$ in the equation for $x_1$:
$$ 2 x_1 - 2 = 0 \implies x_1 = 1. $$
Now we use $x_0,x_1$ in the equation for $x_2$:
$$ 2 x_2 + 1 - 2 = 0 \implies x_2 = \frac{1}{2}. $$
So our solution is
$$ x = x_0 + \epsilon x_1 + \epsilon^2 x_2 + \cdots = 1 + \epsilon + \frac{1}{2}\epsilon^2 + \cdots. $$
If you compare this to the Taylor series expansion of the exact solution $x=\epsilon+\sqrt{1+\epsilon^2}$ you'll see they are exactly the same. :)
In this simple problem perturbation theory is definitely overkill, but in a complicated physical situation it is often essential, and usually the lowest couple of orders are enough accuracy in practice if perturbation theory is going to be useful at all. There are many ways that perturbation theory can go bad, or at least become very tricky, but it is a powerful tool to have in your belt anyway.
