This is potentially a hard problem as it consists of two integral equations. There are also many possible solutions, so you need to break the degeneracy somehow. We can try by using the simplest possible a that can satisfy your general constraints: linearly varying acceleration. In practice I suspect this may not work well at all, but we can try.
Code:
v(t) = v0 + Integral(0,t) a(t') dt'
x(t) = x0 + Integral(0,t) v(t') dt'
v(tf) = v1
x(tf) = x1
a = a0 + a1 t
v(t) = v0 + a0 t + a1 t^2/2
x(t) = x0 + v0 t + a0 t^2/2 + a1 t^3/6
v1 - v0 = a0 t + a1 t^2/2
x1 - x0 - v0 t = a0 t^2/2 + a1 t^3/6
[ x1 - v0 ] = [ t t^2/2 ] [ a0 ]
[ x1 - x0 - v0 t ] [ t^2/2 t^3/6 ] [ a1 ]
That's a system of linear equations you can solve to find one such non constant acceleration. It sucks though, the determinant of that matrix is -t^4/12 which will be a very small number, so this could have serious numerical stability problems. Your accelerations may be extremely bizarre.
When you discover in what way it's bizarre, we can pick an interpolation with more degrees of freedom and a constraint that eliminates the bizarreness.
JASS: