Solving for a cubic polynomial such that the maximum value is constrained over an interval is not simply done, although I do it (approximately) in my SLM tool box.
The trick is to understand the difference between necessary and sufficient constraints on your goal.
Suppose that you picked a finite set of points to test the cubic polynomial at within your interval? If the chosen polynomial satisfies the constraints at each point, then you will accept that this is a necessary condition for the polynomial to be less than your goal everywhere in the interval, recognizing that this is not a sufficient condition. There may still be some points at other locations in the interval where the polynomial exceeds that maximum. Of course your function is only a cubic polynomial, so it must be smooth and continuous. Therefore if you have chosen enough (well placed) points it will not exceed the goal by a lot.
You can define these necessary constraint locations using the linear constraints found in lsqlin, so lsqlin is a viable solver for this. You need not use fmincon at all.
Or, you could just use my slmengine to solve your problem, it defines the set of points as a list of 11 values in the interval, plus it will add the endpoints of the interval to that list. Thus, this call will be all it takes to fit a single cubic polynomial over the interval [0,1], such that it (hopefully) never passes above the value 0.5. It uses the necessary set of constraints I have described above, so it may go slightly above that value though.
x = rand(10,1); y = rand(10,1); slm = slmengine(x,y,'knots',[0 1],'maxvalue',0.5,'result','pp');
This will return a spline in a pp form. Note that the single cubic segment it generates will be relative to the left hand end point, effectively evaluating P(x-x_left).
You can find SLM on the File Exchange:
If you wanted to use more/different point locations, you would need to change the code, but that point is rather easy to find, and easy to modify. (I suppose I could have given that as an option. Ugh.) Or, as I said, just use lsqlin, defining the constraints yourself.
You really cannot easily do any better unless you are willing to get your hands a bit dirty. You would then use fmincon, defining a nonlinear constraint as the maximum value of the cubic over that interval. This is easily found for a cubic. Simply test the points where the cubic polynomial has a zero slope, IF those points lie within the interval of interest. Compare the value found to the value of the cubic at the end points of the interval, returning the overall maximum value as the constraint value. Use of this form of constraint with fmincon will yield a result that is both necessary and sufficient for your goal.
Finally, a flaw with the use of fmincon here is the resulting constraint will not be a differentiable function of the parameters of the polynomial. (The constraint value will be continuous, but it can under some circumstances be non-differentiable.) fmincon will probably survive it, although it may leave you at a sub-optimal location.
