algorithm for convex $C^2$ interpolation

If $c_i$'s are all positive, there are infinitely many such convex $C^2$ functions. As I have pointed out in my comment above, nonnegativity of $c_i$'s is insufficient to guarantee the existence of a $C^2$ function. One very simple construction via Bezier curve is as follows.

Draw a straight line through each point $(x_i,f_i)$ such that all other points lie above the line. In each interval, construct a quartic Bezier curve as follows. Set the control points of the quadratic Bezier curve resulting from the previously drawn straight lines. Then make the midpoint of each line segment a control point (doubling the control number minus one). Draw the quartic Bezier curve from these control points.

The reason for the construction is that the tangent vector (first derivative) with respect to the parameter ($t$ in the Wikipedia article) of a Bezier at an end point is the attached line segment while the second derivative of the curve is the difference of the two closest line segments. We are making the tangent vectors on both sides of the end (data) point coincide and the difference vectors vanish thus equal. You get a $C^2$ curve with a continuous first derivative and a continuous second derivative vanishing at each data point.

The above algorithm, proving the existence of the convex $C^2$ interpolation function, has a vanishing second derivative at each data point. That makes the first derivative run parallel to the $x$ axis each time the curve reaches a data point, making the first derivative wiggly. It does not have to Having proved the existence of an $C^2$ interpolation, we can make the first derivative of the convex $C^2$ interpolating curve smoother by constructing a higher order Bezier curve by connecting neighbouring data points with many small line segments of almost equal lengths each turning almost a constant angle. This will eliminate the horizontal running points from the first derivative and makes it appear smoother.