I'm a bit confused now. My cubic curves have turned out identical to my bezier curves.

I've made a small test program which allows me to specify 4 points to create a path. For the bezier curve, these 4 points [(x0, y0), ..., (x3, y3)] are treated as control points and simply used as inputs to the bezier function. Excuse the slightly messy code...

	private float bezierX(float t){
		float omt = 1 - t; //One minus T
		
		float x01 = mix(x0, x1, t, omt);
		float x12 = mix(x1, x2, t, omt);
		float x23 = mix(x2, x3, t, omt);

		float x012 = mix(x01, x12, t, omt);
		float x123 = mix(x12, x23, t, omt);

		float x0123 = mix(x012, x123, t, omt);
		
		return x0123;
	}
	
	private float bezierY(float t){
		float omt = 1 - t; //One minus T
		
		float y01 = mix(y0, y1, t, omt);
		float y12 = mix(y1, y2, t, omt);
		float y23 = mix(y2, y3, t, omt);

		float y012 = mix(y01, y12, t, omt);
		float y123 = mix(y12, y23, t, omt);

		float y0123 = mix(y012, y123, t, omt);
		
		return y0123;
	}

Now, a bezier curve has the following derivation:

	private float bezierDX(float t){
		float omt = 1 - t; //One minus T
		
		return 3*omt*omt*(x1 - x0) + 6*omt*t*(x2-x1) + 3*t*t*(x3-x2);
	}
	
	private float bezierDY(float t){
		float omt = 1 - t; //One minus T
		
		return 3*omt*omt*(y1 - y0) + 6*omt*t*(y2-y1) + 3*t*t*(y3-y2);
	}

Therefore the initial velocity (velocity at t = 0) is obviously

(3*(x1 - x0), 3*(y1 - y0))

and the final velocity (velocity at t = 1) is

(3*(x3 - x2), 3*(y3 - y2))

Then we can create a cubic parametric curve with the same initial and final velocity (note the values for f3 and f4):

	private float cubicX(float t){
		
		float t2 = t*t;
		float t3 = t2*t;
		
		float f1 = x0 * (1 - 3*t2 + 2*t3);
		float f2 = x3 * (3*t2 - 2*t3);
		float f3 = (x1-x0)*3 * (t - 2*t2 + t3);
		float f4 = (x3-x2)*3 * (-t2+t3);
		
		return f1 + f2 + f3 + f4;
	}
	
	private float cubicY(float t){
		
		float t2 = t*t;
		float t3 = t2*t;
		
		float f1 = y0 * (1 - 3*t2 + 2*t3);
		float f2 = y3 * (3*t2 - 2*t3);
		float f3 = (y1-y0)*3 * (t - 2*t2 + t3);
		float f4 = (y3-y2)*3 * (-t2+t3);
		
		return f1 + f2 + f3 + f4;
	}

Now, the twist here is that these two curves are 100% identical! Assuming I've done everything right so far, I guess I'll just go with cubic parametric curves since it's easier to specify the initial and final velocity.

The problem then is that these cubic curves seem to have widely varying acceleration, just like bezier curves have (since they're identical). This is very problematic since my objects have a maximum acceleration. My guess here is that this can be solved by changing the time given to follow the curve. So by default T goes from 0 to 1, so we'll say that this is in seconds. An object that's going to move 1000 meters in 1 second is going to have to accelerate insanely fast at the beginning to get there in time and then brake extremely fast as well so that it ends up with the correct velocity. If I instead give function 2 seconds to get there, (assuming I've got this right) we need to multiply the initial velocity by 2 since time is passing half as quickly now, but the acceleration becomes 1/4th of the original value.

Does this all make sense?

Perhaps it would be obtain the maxiumum of the acceleration curve by differentation, call it T and divide the curve by (maxAcceleration / T), so that the maximum acceleration will match the object's.

Now I'm not sure, but because integrals are linear you could then simply stretch out the t parameter like you said by the inverse (T/maxAcceleration).

Yes, the curves I'm getting are starting to make sense. There are many interesting scenarios where the algorithm produces some extremely "smart" paths, like where it loops around to be able to gather enough speed without having to accelerate as quickly.

I believe my main problem then is to handle this when under the influence of external forces. Are there any good approaches here?

As Ferrous predicted, I am indeed trying to plan motion for spaceships. External forces would pretty much exclusively be gravity wells (planets, stars). For now all gravity wells are stationary. These are predictable as they only depend on the current position of the object, but to plan/predict the motion of an object I'd most likely need to integrate (e.g. simulate it exactly as I would when actually moving the object).

My biggest problem so far when solving similar problems by iteratively attempting to improve the results is that it's hard to measure the degree of success certain parameters produced.

- Do I terminate the simulation after a set time and measure how off I am from the target position and velocity? I don't know the optimal time.

- Do I terminate the simulation when I get close enough to the target? What if it was going to loop around and get a much better result then?