The given equation is rather cryptic (e.g., $o_i$ are not clearly defined) and thus I will rather address the original problem of finding two palindromes of $d=2l+1$ digits each in bases $b\geq 2$ and $b+1$. This corresponds to solving the equation: $$\sum_{i=0}^{l-1} a_i (b^i + b^{2l-i}) + a_l b^l = \sum_{i=0}^{l-1} c_i ((b-1)^i + (b-1)^{2l-i}) + c_l (b-1)^l$$ in integers $a_0\in[1,b-1]$, $c_0\in[1,b-2]$, $a_i\in [0,b-1]$ and $c_i\in[0,b-2]$ for $i\in\{1,2,\dots,l\}$.
I will show how to solve this equation in a finite number of steps (in particular, finding all finite and infinite series of solutions). For the sake of exposition, let us consider a particular value of $d=5$ ($l=2$).
Step 1. We represent the equation in the form $P=0$, where $P$ is a polynomial in $b$ with coefficients being linear functions in $a_i,c_i$: $$P := (a_0 - 2c_0 + 2c_1 - c_2) + (a_1 + 4c_0 - 4c_1 + 2c_2)b + (a_2 - 6c_0 + 3c_1 - c_2)b^2 + (a_1 + 4c_0 - c_1)b^3 + (a_0 - c_0)b^4.$$
Step 2. We linearize the equation $P=0$ as follows. First, from the bounds for $a_i,c_i$ we obtain bounds for the free term of $P$ (i.e., the coefficient of $b^0$): $$a_0 - 2c_0 + 2c_1 - c_2 \in [1,b-1] - 2[b-2,1] + 2[0,b-2] - [b-2,0] = [-3b+7,3b-7].$$ Then we notice that $P=0$ implies that the free term of $P$ is divisible by $b$, that is $$a_0 - 2c_0 + 2c_1 - c_2 = k_0 b$$ for some integer $k_0$. From the bounds above we have $-3 + \tfrac{7}{b} \leq k_0 \leq 3-\tfrac{7}{b}$, implying that $k_0\in [-2,2]$.
Next, we replace the free term in $P$ with $k_0 b$ and divide the equation $P=0$ by $b$, obtaining $$k_0 + a_1 + 4c_0 - 4c_1 + 2c_2 + (a_2 - 6c_0 + 3c_1 - c_2)b + (a_1 + 4c_0 - c_1)b^2 + (a_0 - c_0)b^3=0.$$ Here we again consider the free term that must be divisible by $b$ and replace it with $k_1b$, and so on.
This results in the system of equations: $$\begin{cases} a_0 - 2c_0 + 2c_1 - c_2 = k_0 b, \\ k_0 + a_1 + 4c_0 - 4c_1 + 2c_2 = k_1b,\\ k_1 + a_2 - 6c_0 + 3c_1 - c_2 = k_2 b,\\ k_2 + a_1 + 4c_0 - c_1 = k_3b,\\ k_3 + a_0 - c_0 = 0, \end{cases} $$ where $k_0\in [-2,2]$, $k_1\in [-3, 6]$, $k_2\in [-6, 3]$, $k_3\in [-1, 4]$.
Step 3. We iterate the $k_i$ over their ranges to obtains a finite number of systems of linear equations over variables $a_i$, $c_i$, and $b$. Together with the bounding conditions for $a_i$ and $c_i$, each such system defines a polyhedron (possibly unbounded), whose integer points can be found with existing algorithms.
For example, this can be done in SageMath with integral_points_generators() function, which uses the PyNormaliz backend.
I implemented this the described algorithm in SageMath, and confirm that the solutions for $d=5$ listed in the table are complete modulo the following typos:
- In the rows labeled $(1,4,4,2)$, the base should be $2a_1+a_2+1$ rather than $a_1+a_2+1$;
- In the last five rows, the value of $a_2$ should be decreased by $1$ (e.g., $2(x+6)$ instead of $2(x+6)+1$).
This way we can get all solutions for $d=7$ and possibly larger $d$'s, but Step 3 needs to be optimized to avoid choices of $k_i$'s that are not feasible. (I'm working on implementing this now.)