If you want to generate the formulas then it is just a matter of solving a set of simultaneous equations
p1=p0+v0*t01+(a0/2)*t01^2+(j/6)*t01^3
v1=v0+a0*t01+(j/2)*t01^2
a1=a0+j*t01
p2=p1+v1*t12+(a/2)*t12^2+(j/6)*t12^3
v2=v1+a*t12+(j/2)*t12^2
a2=a+j*t12
etc until you get to point 7.
a1=a2=a is the max acceleration and j is the max jerk
I think there are 17 or 19 instead of the obvious 21 simultaneous equations to solve because not all segments have a non zero jerk and acceleration. You can do it by hand using substitution.
t01 is the time it takes to go from point 0 to point 1.
Solving for the simultaneous equations is easy. The problem is that there are so many cases. For instance the motion profile may not reach the max acceleration or deceleration. It may not even reach the maximum velocity in the case of short moves. Then what. What happens when you issue a new command while ramping up or ramping down or the position is now in the opposite direction from the current one.
After all the positions, velocities, accelerations and times are work out you can build a state machine with a switch or case statement where each state has the position, velocity and acceleration formulas for each segment.