foresightful acceleration and decelerations

I am programming a CNC machine (using matlab). In order to generate a surface with a "high" shape accuracy and micro roughness I need to make sure that no large acceleration/decelerations takes place, because this would induce vibrations. Currently I use a moving average, which I apply to the distance of the data-point themselves. Although the result seams promising, I don’t like this approach conceptionally. Instead I would love to apply the moving average to the acceleration itself and then use $ $ x(i+1) = x(i) + v(i) dt + \frac{1}{2} \bar{a}(i) dt^2 $ $ However, my problem with this approach is that the code gets messy, because I have several axes, each has a maximal velocity and maximal acceleration. Are there other simple methods available? How would you implement the above?