How to create State-space equation from complex symbolic equation?

I generated a code for robotic dynamics with lagrange euler method. I generated for 4 dof. I want to set up the torque equations with state-space matrix definitions.
I found the thau torque for one joint like;
Izz1*dd_the1 + Iyy2*dd_the1*cos(the2)^2 + Iyy4*dd_the1*cos(the2)^2 + Izz3*dd_the1*cos(the2)^2 +
Ixx2*dd_the1*sin(the2)^2 + Iyy3*dd_the1*sin(the2)^2 + Izz4*dd_the1*sin(the2)^2 - (dd_the1*l2^2*m2
*cos(the2)^2)/4 - (dd_the1*l4^2*m4*cos(the2)^2)/4 + (d_the2^2*l3^2*m3*sin(the2))/4 + (d_the2^2*l3
^2*m4*sin(the2))/2 - (dd_the1*l2^2*m2*sin(the2)^2)/4 - (dd_the1*l3^2*m3*sin(the2)^2)/4 - (dd_d4
*l3*m4*cos(the2))/2 - (dd_d3*l3*m3*sin(the2))/2 - (dd_d3*l3*m4*sin(the2))/2 + (d_the1^2*l3^2*m3
*sin(the2)^3)/4 + (d_the1^2*l3^2*m4*sin(the2)^3)/2 - (dd_the2*l3^2*m3*cos(the2))/4 - (dd_the2*l3^
2*m4*cos(the2))/2 + 2*Ixx2*d_the1*d_the2*cos(the2)*sin(the2) - 2*Iyy2*d_the1*d_the2
*cos(the2)*sin(the2) + 2*Iyy3*d_the1*d_the2*cos(the2)*sin(the2) - 2*Iyy4*d_the1*d_the2
*cos(the2)*sin(the2) - 2*Izz3*d_the1*d_the2*cos(the2)*sin(the2) + 2*Izz4*d_the1*d_the2
*cos(the2)*sin(the2) + (d4*d_the2^2*l3*m4*cos(the2))/2 - (d4*dd_the1*l4*m4*cos(the2)^2)/2 - d_d3
*d_the1*l4*m4*cos(the2)^2 - (d_the1^2*l2*l4*m4*cos(the2))/2 + (d_the2^2*l3*l4*m4*cos(the2))/4 -
(dd_the1*l2*l3*m3*cos(the2)^2)/2 - (dd_the1*l2*l3*m4*cos(the2)^2)/2 + (d3*d_the2^2*l3*m3
*sin(the2))/2 + (d3*d_the2^2*l3*m4*sin(the2))/2 - (d3*dd_the1*l3*m3*sin(the2)^2)/2 - (d3*dd_the1
*l4*m4*sin(the2)^2)/2 - d_d4*d_the1*l4*m4*sin(the2)^2 - (g*l4*m4*cos(the1)*cos(the2))/2 - (d_the1
^2*l2*l3*m3*sin(the2))/2 - (d_the1^2*l2*l4*m4*sin(the2))/2 - (dd_the1*l2*l3*m3*sin(the2)^2)/2 -
(dd_the1*l2*l3*m4*sin(the2)^2)/2 - (dd_the1*l3*l4*m4*sin(the2)^2)/2 - (g*l3*m3
*cos(the1)*sin(the2))/2 - (g*l4*m4*cos(the1)*sin(the2))/2 + (d_the1^2*l3^2*m3*cos(the2)^2
*sin(the2))/4 + (d_the1^2*l3^2*m4*cos(the2)^2*sin(the2))/2 + (d4*d_the1^2*l3*m4*cos(the2)^3)/2 +
(d_the1^2*l3*l4*m4*cos(the2)^3)/4 + (d3*d_the1^2*l3*m3*sin(the2)^3)/2 + (d3*d_the1^2*l3*m4
*sin(the2)^3)/2 + (d_the1*d_the2*l4^2*m4*sin(the2)^2)/2 - (dd_the1*l4^2*m4
*cos(the2)*sin(the2))/4 - (g*l2*m2*cos(the2)^2*sin(the1))/2 - (g*l3*m3*cos(the2)^2*sin(the1))/2 -
(g*l3*m4*cos(the2)^2*sin(the1))/2 - (g*l2*m2*sin(the1)*sin(the2)^2)/2 - (g*l3*m3
*sin(the1)*sin(the2)^2)/2 - (g*l3*m4*sin(the1)*sin(the2)^2)/2 - (d3*dd_the2*l3*m3*cos(the2))/2 - (d3
*dd_the2*l3*m4*cos(the2))/2 + (d4*dd_the2*l3*m4*sin(the2))/2 + d_d4*d_the2*l3*m4*sin(the2) +
(dd_the2*l3*l4*m4*sin(the2))/4 - (d_the1*d_the2*l3^2*m3*cos(the2)*sin(the2))/2 + (d_the1*d_the2
*l4^2*m4*cos(the2)*sin(the2))/2 + (d3*d_the1^2*l3*m3*cos(the2)^2*sin(the2))/2 + (d3*d_the1^2*l3
*m4*cos(the2)^2*sin(the2))/2 + (d4*d_the1^2*l3*m4*cos(the2)*sin(the2)^2)/2 + (d_the1^2*l3*l4*m4
*cos(the2)*sin(the2)^2)/4 - d3*d_the1*d_the2*l4*m4*cos(the2)^2 - d_the1*d_the2*l3*l4*m4
*cos(the2)^2 + d4*d_the1*d_the2*l4*m4*sin(the2)^2 - (d3*dd_the1*l4*m4*cos(the2)*sin(the2))/2 -
(d4*dd_the1*l4*m4*cos(the2)*sin(the2))/2 + d_d3*d_the1*l3*m4*cos(the2)*sin(the2) - d_d3*d_the1
*l4*m4*cos(the2)*sin(the2) - d_d4*d_the1*l4*m4*cos(the2)*sin(the2) - (dd_the1*l3*l4*m4
*cos(the2)*sin(the2))/2 - d3*d_the1*d_the2*l3*m3*cos(the2)*sin(the2) - d3*d_the1*d_the2*l4*m4
*cos(the2)*sin(the2) + d4*d_the1*d_the2*l4*m4*cos(the2)*sin(the2) - d_the1*d_the2*l3*l4*m4
*cos(the2)*sin(the2)

How can i create state-space from this. Or Is there a matlab algorithm that looks for coefficients of any specific variable from the symbolic equation

The state-space form i want to create;

https://i.hizliresim.com/7r47Ye.png

[delil codes]

Variables;

dd_the , dd_d is for D(q)*qDotDot

g is for G(q)

d_the , d_d is for C(q,qDot)