This is a very nice question - an innocent-looking problem, easily described - basically construct a horizontal s-shaped analytical function under some well-defined conditions.
But I must admit I had changed my mind several times, bouncing back and forth, in trying to decide (and prove) whether it can or cannot be done.
I now believe that the problem can be solved with polynomials as @MattF wrote in his answer,
and as you and @LoïcTeyssier suggested in the comments, the degree $n$ of the polynomial depends on the constraints.
The problem is very suited to be represented with polynomials in Bernstein/Bezier form (see also this comprehensive survey).
The Bernstein polynomial $P(x)$ is defined as $P(x) = \sum_{i=0}^{n} c_i B_{i}^{n}(x)$,
where $B_i^{n}$ are the Bernstein basis polynomials $B_i^{n}(x) = \binom{n}{i} (1-x)^{n-i} x^{i}$.
Below, I will sometimes use freely the term "curve" (or "polynomial curve" or "Bezier curve") for the graph of the polynomial.
Also, the points $p_i = (i/n, c_i)$, denoted the control points or control polygon of the curve, have a geometric meaning, which I will be using.
The first thing to notice is that the conditions on the first and second derivatives totally define the first and last three control points.
Since the Bezier curves interpolate their endpoints, we have from condition (0):
$$c_0 = c_n = 0$$
The first derivative of a Bernstein polynomial of degree $n$ is itself a Bernstein polynomial of degree $n-1$ of the form:
$P'(x) = n \sum_{i=0}^{n-1} (c_{i+1} - c_i) B_{i}^{n-1}(x)$
Thus $P'(0) = n (c_1 - c_0)$ and from condition (2) we get:
$n (c_1 - c_0) = f'(0)$,
and so:
$$c_1 = f'(0)/n$$
and in a similar manner:
$$c_{n-1} = -f'(1)/n$$.
The second derivative is a Bernstein polynomial of degree $n-2$ of the form:
$$P''(x) = n(n-1) \sum_{i=0}^{n-2} (c_{i+2} - 2 c_{i+1} + c_i) B_{i}^{n-2}(x)$$
Therefore,
$$P''(0) = n(n-1)(c_2 - 2 c_1 + c_0)$$
and from condition (2) we get:
$$n(n-1) (c_2 - 2 c_1) = f''(0)$$
and (with some manipulation):
$$c_2 = \frac{f''(0)}{n(n-1)} + 2 c_1 = \frac{f''(0)}{n(n-1)} + 2 \frac{f'(0)}{n}$$.
and in a similar manner:
$$c_{n-2} = \frac{f''(1)}{n(n-1)} - 2 c_{n-1} = \frac{f''(1)}{n(n-1)} - 2 \frac{f'(1)}{n}$$.
Furthermore, if the control polygon intersects the $x$-axis only once in $(0, 1)$ then, by the Bezier variation diminishing property, there is only one inner root.
This is useful since we want (from condition (0)) only a single inner root (at $x=p$),
and the roots at $0$ and $1$ are immediate since we set $c_0 = c_n = 0$.
Another important concept is what I call the "convex-concave property" of the control polygon.
I say the control polygon is "convex-concave" if every consecutive control point
tuple $(p_i, p_{i+1}, p_{i+2})$ make a left-turn up to some index $i=k$ and from then on all consecutive tuples make right-turns.
The important insight here is that if the control polygon is convex-concave, then so is its Bezier curve.
The reason is that, as we saw above, the second derivative of the polynomial is itself a Bezier curve (of degree $n-2$), where its coefficients are $d_i = c_{i+2} - 2c_{i+1} + c_i$.
A left-turn tuple gives a positive control point $d_i$, whereas a right-turn tuple gives a negative one.
Thus, up to $d_k$ all $d_i$ are positive, and afterwards all are negative.
Therefore, by the variation diminishing property,
there is only one zero-crossing of the second-derivative control polygon,
which means that there is just a single inflection point in $[0, 1]$.
So a convex-concave control polygon certifies that the curve is convex up to some point and concave afterwards.
One thing to notice is that if $c_2 > c_1$ then it may be impossible to satisfy the convex-concave condition.
To see this, take for example $c_2 = 0$. Then $c_3$ must be positive or the tuple $(c_1, c_2, c_3)$ will be concave and violate the convex-concave property.
If this happens on both sides (e.g., $c_2 = 0$ and $c_{n-2} = 0$) then the convex-concave property cannot be satisfied.
If, however, $c_2 < c_1$ and $c_{n-2} > c_{n-1}$, then you can always construct a convex-concave polygon.
For example, we can construct such a polygon by taking all $c_i$ to be on
the line between $c_1$ and $c_2$ for $i < k$,
and all $c_i$ on the line between $c_{n-1}$ and $c_{n-2}$ for $i \geq k$.
We can derive a lower bound on $n$ for which this will be the case.
Since $c_2 = \frac{f''(0)}{n(n-1)} + 2 c_1 = \frac{f''(0)}{n(n-1)} + 2 \frac{f'(0)}{n}$,
then if we set $n > |f''(0)/f'(0)| + 1$, we will have $c_2 < c_1$.
Similarly, for $n > |f''(1)/f'(1)| + 1$, we will have $c_{n-2} > c_{n-1}$
So setting the lower bound $n > \max(|f''(0)/f'(0)| + 1, |f''(1)/f'(1)| + 1)$
guarantees that we can construct convex-concave control polygons with the required end conditions.
This already takes us a long way.
Condition (2) is satisfied by the first and last three control points.
Condition (1) is satisfied in the weak sense (i.e., with an implicit $p$) by the convex-concave property.
And condition (0) is also partially satisified by the first and last control points and the fact that there is only one inner root.
What is missing to fully satisfy all the conditions (at least in the weak sense)
is to find a way to force the root and the inflection point
(where the curve turns from convex to concave) to be the same point $p$.
We can now think of an algorithm to construct the polynomial, which starts with the first and last three control points,
and continues by trying to build a convex-concave control polygon, which satisfies the additional constraints.
Apart from the constraints on the first and last three coefficients described above, there are the following constraints.
The polynomial zero-value at $p$:
$$P(p) = \sum_{i=0}^{n} c_i B_{i}^{n}(p) = 0$$
The zero-value of the second derivative at $p$:
$$P''(p) = \sum_{i=0}^{n-2} (c_{i+2} - 2 c_{i+1} + c_i) B_{i}^{n-2}(p) = 0$$
The inequality constraints to satisfy the convex-concave property:
$$(c_{i+2} - 2 c_{i+1} + c_i) > 0 \text{, for } 0 < i < k$$
$$(c_{i+2} - 2 c_{i+1} + c_i) < 0 \text{, for } k \leq i < n-2$$
These constraints are linear equalities and inequalities in the variables $c_i$ and can therefore be solved using a Linear Programming solver (with varying k-values).
Note that the above lower bound on $n$ may not be sufficient.
For example, if $p$ is very close to $0$ or $1$, e.g., smaller than $1/n$, then
we will need to increase $n$.
However, the LP-algorithm also reports if a feasible solution is not found and one can then increase $n$ and try again.
Even though I believe this algorithm is theoretically correct (after all, there are so many degrees of freedom..) I haven't succeeded in proving it.
I wrote an initial implementation of the algorithm, and while it works fairly well on some inputs, on others the LP solver fails due to numerical errors (I used scipy's LP solver).
This looks to me like an implementation issue (possibly a bug in my code), but there might be something fundamental hiding there, which I'm missing.
Below is a figure of polynomials and their second derivatives for $p=0.5$ and varying end conditions, computed by the algorithm.
For the weaker version of the problem (i.e., with an implicit $p$), I believe I can show that a solution exists.
Given the first and second derivative constraints on the endpoints, I'll show that there exists a convex-concave polynomial $P(x)$ (of a sufficiently large degree $n$) that satisfies the end constraints,
and for which there exists some value $p$ that satisfies $P(p)=0$ and $P''(p)=0$.
This means that $x=p$ is both the root and the inflection point and with the convex-concave property satisfies the
weaker conditions of the problem. Note, I prove the existence of $p$ for a sufficiently large $n$, but don't explicitly compute $p$ nor $n$.
I will use the Weierstrass approximation theorem, which states that any continuous function on a closed interval can be approximated with a polynomial.
More specifically, I will use Bernstein's constructive proof of the Wierstrass theorem, which gives the following theorem:
For a continuous function $f$ over $[0, 1]$, for any $\epsilon$, for a sufficiently large $n$,
$|f(x) - \sum_{i=0}^n f(i/n) B_{i}^{n}(x)| < \epsilon$.
From Bernstein's theorem, by setting the control polygon to be the function $f$, the following corollary follows:
For any $\epsilon$, for a sufficiently large $n$, the distance between the Bernstein polynomial and its control polygon is smaller than $\epsilon$.
From the corollary it follows that (for a suitable $n$)
the root $p$ of $P(x)$ is within a $\delta$-environment of the intersection point of the control polygon and the $x$-axis.
Similarly, the inflection point (the root of $P''(x)$) is within a $\delta$-environment of the intersection of the control polygon of $P''(x)$ and the $x$-axis.
For any $n$, I assume we can construct a convex-concave control polygon, which satisfies the end conditions,
and has an inflection point above the $x$-axis, i.e., to the right of the root.
For a sufficiently large $n$, by the corollary, the curve inflection point will also be to the right of the curve root.
We can similarly construct such a configuration that has an inflection point to the left of the root.
Below is are example figures of two such corresponding configurations, where the inflection point is the blue dot, and the root is black (the control points are the blue crosses).
Moving continuously from the first configuration to the second (while maintaining the convex-concave property),
we will encounter (at least one) value $x=p$ where the inflection
point and the root coincide, which satisfies our request.
So, to summarize, I assume we can draw for any $n$ a pair of convex-concave control polygons that satisfy the end conditions and for which one has an inflection point above the $x$-axis and the other has it below (I think this is a reasonable assumption).
Then, for a sufficiently large $n$, the polynomial curves will be close enough to their control polygon so that their inflection point is also to the right/left of the polynomial root.
Then there should be at least one configuration where the root and inflection point coincide.
Some open issues:
Try and convert the weak-constraint proof to an algorithm. Basically, if you have the two configurations with their roots on opposite sides of the inflection point, then a bisection algorithm might be used.
Prove that there is always a solution for the problem with the strong constraints (for a given $p$ value).
The consequence is that the LP algorithm will eventually give the correct solution.
Find a (hopefully analytical/concise/easily computable) upper bound on the polynomial degree $n$.
