Form sub-equations from a parent equation

I have 2 coupled differential equations which I want to solve using Galerkin Method. As such, I want to substitute the dependent variables as a polynomial expansion of the independent variable, meantime satisfying the boundary conditions. Let me put the thing mathematically-

The 2 differential equations are-
$$
frac{partial x}{partial t}=f_1(x,y,t)
frac{partial y}{partial t}=f_2(x,y,t)
$$
To solve it using Galerkin Method, I use the following polynomial expansion-
$$
x(t) = a_1t(1-t) + a_2t^2(1-t) + … = Sigma_{i=1}^n a_it^i(1-t)
y(t) = b_1(1-t) + b_2(1-t^2) + … = Sigma_{i=1}^n b_i(1-t^i)
$$
Substituting the polynomials in the governing differential equations yields a system of coupled algebraic equations with the coefficients $a_i$ and $b_i$ forming the unknown variables.

My question is how to automate the process of forming the algebraic equations after substituting the polynomials in the differential equations? Use suitable functions to simplify the problem. I just need to know how to find the algebraic equations.