You can use a fifth-degree (quintic) polynomial to describe motions for each axis. It is also possible to use other polynomials (second-degree, third-degree, etc.) but for polynomials with degrees lower than five, you will need a concatenation of polynomials (i.e. piecewise-polynomials) to describe motions of each axis.
Suppose the motions you are interested in starts and ends at rest (i.e. having zero velocity and acceleration). With this condition, you only need to calculate only one set of quintic coefficients $a, b, c, d, e, f$. Then the quintic cofficients for each of the axis will be just a scaled version of those coefficients.
To do this, let's suppose we have an extra, imaginary axis for this robot (just for this calculation) that moves from $0$ to $1$ and whose velocity limit $v_m$, acceleration limit $a_m$, and maybe jerk limit $j_m$ are determined by
$$
\begin{align}
v_m &= \min_{i} v_{m,i}/\text{displacement}_{i},\\
a_m &= \min_{i} a_{m,i}/\text{displacement}_{i},\\
j_m &= \min_{i} j_{m,i}/\text{displacement}_{i},
\end{align}
$$
where $\text{displacement}_{i}$ is how much axis $i$ moves, $v_{m,i}$ is the max velocity of axis $i$, and so on.
Now we are going to determine an appropriate quintic polynomial for this axis:
$$
p(t) = at^5 + bt^4 + ct^3 + dt^2 + et + f,
$$
which will describe a motion from time $t = 0$ to $t = T$, where $T$ is the total motion duration, which is still not determined at this point.
After we finish the calculation, we can then compute quintic coefficients for all "actual" robot axes just by scaling the above $a, b, \ldots, f$.
Given the boundary conditions (zero velocity and acceleration), we can compute the coefficients as follows
$$
a = \frac{6}{T^5}, b = \frac{-15}{T^4}, c = \frac{10}{T^3}, d = 0, e = 0, f = 0.
$$
So now we can actually view the polynomial equation as parameterized by two variables $p = p(t, T)$.
If there were no velocity/acceleration limits, etc. for the motor, we can assign $T$ arbitrarily. However, with the limits, we have to make sure that the computed polynomial is respecting the limit.
I'll just show here how to select $T$ based on the velocity limit $v_m$. Otherwise, it's going to be too long.
From the polynomial equation, the velocity profile can be computed as
$$
v(t, T) = \frac{d}{dt}p(t, T).
$$
This velocity profile has its peak at $t_p$, which solves the following equation
$$
\frac{d}{dt}v(t, T) = 0.
$$
Once we solve for $t_p$, we can then get the peak velocity as $v_p = v(t_p, T) = (15/8)/T$. So if we want the peak velocity to reach the maximum, then
$$
v_p = v_m = (15/8)T.
$$
Therefore, if we set $T = (15/8)/v_m$, we will have a quintic polynomial whose velocity profile reaches the max velocity. Call this duration $T_v$. In a similar fashion, we can compute, for example, $T_a$, which is the duration that make the velocity profile reaches max acceleration. Finally, select the duration as $T = \max(T_v, T_a, T_j)$.
Now that we have completely determine the quintic coefficients for this imaginary axis, we can obtain coefficients for axis $i$ from $a_i = a*\text{displacement}_i$, $b_i = b*\text{displacement}_i$, and so on.
