Question: Let $f$ be a real valued function defined on an interval $I$. Prove that if for any $x,y,z\in I$ satisfying $x<y<z$ we have $$\begin{vmatrix} 1 & 1 & 1\\ x & y & z\\ f(x) & f(y) & f(z)\end{vmatrix}\ge 0,$$ then $f$ is convex.
Solution: Observe that for any $x\neq y\in I$, $x=\lambda x+(1-\lambda) y\iff \lambda =1$ and $y=\lambda x+(1-\lambda)y\iff \lambda=0$.
Now let us take any $x\neq y$ such that WLOG we have $x>y$. Also select any $\lambda\in(0,1)$.
Thus, we have $$\lambda x> \lambda y\implies \lambda x+(1-\lambda) y>\lambda y+(1-\lambda)y=y.$$
Again since $\lambda<1\implies -\lambda>-1\implies 1-\lambda>0.$ Therefore since $$x>y\implies (1-\lambda)x>(1-\lambda)y\implies x=\lambda x+(1-\lambda)x>\lambda x+(1-\lambda)y.$$
Thus, for any $\lambda\in(0,1)$ we have $$x>\lambda x+(1-\lambda)y>y.$$
This implies that by our hypothesis we must have $$D:=\begin{vmatrix}1 & 1 & 1\\ x & \lambda x+(1-\lambda)y & y\\ f(x) & f(\lambda x+(1-\lambda)y) & f(y) \end{vmatrix}\ge 0.$$
Applying the column transformations $C_2\to C_2-\lambda C_1$, followed by $C_2\to C_2-(1-\lambda)C_3$, we have $$D=\begin{vmatrix}1 & 0 & 1\\ x & 0 & y\\ f(x) & f(\lambda x+(1-\lambda)y)-\lambda f(x)-(1-\lambda)f(y) & f(y) \end{vmatrix}\ge 0\\\implies f(\lambda x+(1-\lambda)y)\le \lambda f(x)+(1-\lambda)f(y).$$
Hence, we are done in this case.
Now let $\lambda=0$. Thus $f(\lambda x+(1-\lambda)y)=f(y)=\lambda f(x)+(1-\lambda)f(y)$.
Now let $\lambda=1$. Thus $f(\lambda x+(1-\lambda)y)=f(x)=\lambda f(x)+(1-\lambda)f(y)$.
Finally select any $x,y\in I$ such that $x=y$. Also fix $\lambda \in[0,1]$. Thus, $f(\lambda x+(1-\lambda)y)=f(x)=f(y)=\lambda f(x)+(1-\lambda)f(y).$
Hence, we are done with all the cases and we can conclude that in all cases we have $$f(\lambda x+(1-\lambda) y)\le \lambda f(x)+(1-\lambda)f(y), \forall x,y\in I\text{ and }\forall \lambda\in [0,1].$$
Therefore, we can conclude that $f$ is convex.
Is this solution correct and rigorous enough? And, is there any other way to solve the same?
Here is, as you ask, another way to solve the issue :
$$\dfrac12 \begin{vmatrix} x_1 & x_2 & x_3\\ y_1 & y_2 & y_3\\ 1 & 1 & 1\\ \end{vmatrix}$$
is the oriented area of triangle with vertices $A_k(x_k,y_k) \ \ (k=1,2,3)$, oriented meaning $>0$ is the sequence $A_1A_2A_3$ has a positive (anti-clockwise) orientation, $<0$ otherwise.
Therefore, it suffices to consider the same determinant with
$$y_k=f(x_k) \ \ \text{for all} \ \ k $$
and rotation of rows $1\rightarrow 2 \rightarrow 3 \rightarrow 1$ (which doesn't change the value of the determinant, and in particular its sign) to conclude at the given property (that can even be considered as a strict convexity criteria, if valid for any $x<y<z$).
Reference : among some others, the answer to this question : https://math.stackexchange.com/q/2130167