Using polynomial ease curves, it's easy to create a function f(x) which obtains f(0) = 0 and f(1) = 1 as well as reaching 0 derivatives f'(0) = 0, f'(1) = 0, f"(0) = 0, f"(1) = 0, etc. One can start with a simple cubic polynomial: f(x) = 3x²-2x³ which achieves first order smoothness (i.e. f'(0) = 0 and f'(1) = 0. Then we can compose this to have higher and higher continuity f(f(x)), f(f(f(x))). This creates a 9th and 27th order polynomials. Here's a sequence of repeated composing this cubic function. The title shows the order: We can also derive such a smooth, symmetric polynomial ease curve for any order 2*k+1, k in ℕ^{+} polynomial. Thus with k we can have finer, albeit incremental, control over the shape of the function. Here's a ruby script to generate a maple script that will generate the necessary coefficients (there must be a pattern here, but I'm too lazy right now to find it):

```
#!/opt/local/bin/ruby
# Script to generate Maple code for solving for 2*k+1 order polynomial f(x)
# with value f(0) = 0, f(1) = 2 and d<sup>k</sup>f/dx<sup>k</sup>|<sub>0</sub> = 0 d<sup>c</sup>f/dx<sup>c</sup>|<sub>1</sub> = 1, c=1,...,k
smoothness = 3
(1..smoothness).each do |k|
as = (k+1..2*k+1).map{|c| "a#{c}"}.join(",");
puts "f := (#{as},x) -> "
puts " #{(k+1..2*k+1).map{|c| "a#{c}*x^#{c}"}.join("+")};"
puts "solve({"
puts " f(#{as},0)=0,"
puts " f(#{as},1)=1,"
diffs = (1..k).map do |c|
" eval(diff(f(#{as},x),#{(1..c).map{|d| "x"}.join(",")}),x=0)=0,\n"+
" eval(diff(f(#{as},x),#{(1..c).map{|d| "x"}.join(",")}),x=1)=0"
end
puts "#{diffs.join(",\n")}},"
puts " {#{(k+1..2*k+1).map{|c| "a#{c}"}.join(",")}});"
end
```

Using these coefficients we can plot the first few: These sequences converge on a scaled and shifted Heaviside step function. But we only have an integer-valued parameter. It'd be nice to have a smooth parameter. If we only cared about this property of the curve getting steeper in the middle as we increase or parameter (while maintaining interpolation and at least C^{1} continuity, then we could try to do this with a cubic spline. Now we can smoothly increase the steepness, but before reaching the step function the spline curves to far, losing injectivity when treated as a function of x. The fact the we cannot reproduce the Heaviside function with a cubic spline should have been obvious anyway. The Heaviside function is *loosely* a degree ∞ polynomial, with ∞ smoothness at 0 and 1. We can achieve what we want by taking a function which indeed converges to the Heaviside step function: f(x) = 1/(1+e^{-2tx}). We can adapt this to our unit square domain and enforce interpolation for any k: g(x) = (f(2x-1)-1/2)/(f(1)-1/2)/2+1/2, or in MATLAB:

```
g = @(x,t) (1./(1+exp(-2*t*(2*x-1)))-1/2)/(1./(1+exp(-2*t*1))-1/2)/2+1/2);
```

This function behaves pretty nicely for t in (0,∞), going from a linear function to the Heaviside step function: To achieve kth order continuity at our endpoints we can just run this function through (i.e. compose with) one of the 2k+1 order polynomials above. That way, with t we smoothly span a whole space of kth order smooth ease curves. For example, if h(x) = 3x²-2x³ then this is h(g(x)):