Assumptions
- $A \in M_{m\times n}(\Bbb R)$ ($m\le n$) has rank $n$ and basis matrix $B$.
- $x_B$ denotes the basic solution.
- $c_B$ denotes the reduced objective function so that $c^T x=c_B^T x_B$. (So the order/arrangement of basic variables is very important.)
The initial simplex tableau for your problem(, which is already in canoncial form,) will be
\begin{array}{cc|c}
A & I & b \\ \hline
-c^T & 0 & 0
\end{array}
We multiply the $m$ constraints by $B^{-1}$. Note that $\mathbf{e}_1,\dots,\mathbf{e}_n$ can be found in $[B^{-1}A\;B^{-1}]$ since $B$ is chosen from columns of $[A\;I]$. We add a linear combination of the first $m$ row to the last row.
\begin{array}{cc|c}
B^{-1}A & B^{-1} & B^{-1}b \\ \hline
c_B^T B^{-1}A-c^T & c_B^T B^{-1} & c_B^T B^{-1}b \tag{*} \label{t1}
\end{array}
Why is $c_B^T$ chosen? We write $x_B$ as $(x_{B_1},\dots,x_{B_n})^T$. Choose an index $i$ in $\{1,\dots,n\}$, and suppose that $x_j = x_{B_i}$
Case 1: $1\le j \le n$. In the $j$-th column of the last row,
$$(c_B^T B^{-1}A-c^T)_j = c_B^T B^{-1} \mathbf{a}_j -c_j = c_B^T \mathbf{e}_i - c_j = c_{B_i} - c_j = 0.$$
Case 2: $n+1 \le j \le m+n$. Note that $c_{B_i}= 0$. In the $j$-th column of the last row,
$$(c_B^T B^{-1})_{j - n} = c_B^T \mathbf{e}_i = c_{B_i} = 0.$$
Therefore, the entries corresponding to the basic variables in the last row in tableau $\eqref{t1}$ will be zero. Then, $\eqref{t1}$ is a simplex tableau that you can compute from the given optimal solution.
Computational results for your example
Here's the GNU Octave code for finding the optimal tableau. The syntax should be the same in MATLAB except for the solver glpk. You may test the code at Octave Online.
c=[60 100 80]'; A=[144 192 240; 100 150 120; 1 1 1];
b=[120000 60000 500]';
[x_max,z_max]=glpk(c,A,b,[0 0 0]',[],"UUU","CCC",-1)
x=[x_max; b-A*x_max] % Include slack variables
A1=[A eye(3)]; % Part of initial tableau
c1=[c' zeros(1,3)]';
B=A1(:,find(x > 0)')
c_B=c1(find(x > 0));
disp('Initial tableau');
[A1 b; -c1' 0]
disp('Final tableau');
[B^-1*A B^-1 B^-1*b; c_B'*B^-1*A-c' c_B'*B^-1 c_B'*B^-1*b]
The computer output
x_max =
0.00000
400.00000
0.00000
z_max = 4.0000e+04
x =
0.0000e+00
4.0000e+02
0.0000e+00
4.3200e+04
-7.2760e-12
1.0000e+02
B =
192 1 0
150 0 0
1 0 1
Initial tableau
ans =
144 192 240 1 0 0 120000
100 150 120 0 1 0 60000
1 1 1 0 0 1 500
-60 -100 -80 -0 -0 -0 0
Final tableau
ans =
Columns 1 through 6:
6.6667e-01 1.0000e+00 8.0000e-01 0.0000e+00 6.6667e-03 -0.0000e+00
1.6000e+01 0.0000e+00 8.6400e+01 1.0000e+00 -1.2800e+00 0.0000e+00
3.3333e-01 1.1102e-16 2.0000e-01 0.0000e+00 -6.6667e-03 1.0000e+00
6.6667e+00 0.0000e+00 0.0000e+00 0.0000e+00 6.6667e-01 0.0000e+00
Column 7:
4.0000e+02
4.3200e+04
1.0000e+02
4.0000e+04
The optimal tableau:
\begin{equation*}
\begin{array}{ccccccc|l}
& x_1 & x_2 & x_3 & s_1 & s_2 & s_3 & \\ \hline
x_2 & \frac23 & 1 & \frac45 & 0 & \frac{1}{150} & 0 & 400 \\
s_1 & 16 & 0 & \frac{432}{5} & 1 & -\frac{32}{25} & 0 & 43200 \\
s_3 & \frac13 & 0 & \frac15 & 0 & -\frac{1}{150} & 1 & 100 \\ \hline
& \frac{20}{3} & 0 & 0 & 0 & \frac23 & 0 & 40000 \\
\end{array}
\end{equation*}
Note that the "unstable operation" above won't give you a different solution. You've just replaced the basic variable $x_1 = 0$ by a nonbasic variable $x_4 = 0$, but the solution shown in the above optimal simplex tableau is still $(x_1,x_2,x_3,x_4,x_5,x_6)^T = (0,2,3,0,0,0)^T$.
Therefore, only "stable operations" with a non-zero variable $x_r \ne 0$ to be replaced by $x_j$ will give you a different basic optimal feasible solution.
It's easy to verify that any convex combination of a set of basic optimal feasible solution(s) is still an optimal feasible solution (since the feasible region in a linear program in convex), so the set of optimal feasible solution is convex (i.e. path connected). Hence, the answer to your question is yes.
(Added in response to the edit of the question)
You may just think of the graph so as to visualize your problem. A basic solution corresponds to an extreme point in the feasible region. (If $\mathbf x$ is an extreme point, then there doesn't exists two different $\mathbf u$ and $\mathbf v$ such that $\mathbf x$ is a convex combination of $\mathbf u$ and $\mathbf v$.) Form the graph of the set of all optimal feasible solutions of the linear program, and note its convexity (so it's a convex polygon), then eliminate its relative interior point. Take any two vertices of the remaining boundary points to be $X_1$ and $X_2$ in (1). Clearly, we can see a path running from $X_1$ to $X_2$ through adjacent vertices. (Assume that you have $m$ constraints and $n$ decision varibles, and $m < n$. Choose $n$ vectors in the $m$-by-$m+n$ matrix $A$ to form a basis matrix $B$. This is similar to choosing hyperplanes (representing the constraints) in $\Bbb R^{m+n}$ and find its intersection.) This corresponds to a (finite) series of steps of replacing $x_r$ by $x_j$, which is feasible. Since that's feasible, there exists a pivot operation in one of the simplex tableaux $T_1$ for $X_1$.
Let $y_{00} = \bar c$ be the objective function value, $y_{0j}$ be the $j$-th column of the objective function row (i.e. the row of $\bar c$), $y_{r0} = x_r$ be the $r$-th component of the current solution, and $y_{rj}$ be the $r,j$-th entry in the coefficient matrix.
\begin{matrix}
y_{0j}& y_{00}\\
y_{rj}& y_{r0}
\end{matrix}
In order not to change the value of $y_{00}$ after a pivot operation $y_{00}-\dfrac{y_{r0} y_{0j}}{y_{rj}}$, i.e. $$y_{00}-\dfrac{y_{r0} y_{0j}}{y_{rj}} = y_{00},$$ we have $y_{r0} y_{0j} = 0$. Since we want to move to another point, we must have $y_{r0} \ne 0$, which implies $y_{0j} = 0$, so the operation must be stable. Hence (1) is proved.
Best Answer
In my opinion, your constraints do not represent the situation. You can have more than one optimal basic solution if the slope of the objective function is equal to one of the constraints. Then the objective function lies on a line, between two corners of the feasible region. I give you an example with the corresponding graph:
$\textrm{max} \ Z=500x_{1}+300x_{2}$
s.t.
$15x_{1}+5x_{2}\leq 300$
$10x_{1}+6x_{2}\leq 240$
$8x_{1}+12x_{2}\leq 450$
$x_{1},x_{2}\geq 0$
The corresponding graph is
To make the graph we have to solve the (in-)equalities for $x_2$. I focus on the the objective function and on the constraint with the same slope.
1. $\textrm{Objective function}$
$Z=500x_{1}+300x_{2}$
$Z-500x_{1}=300x_{2}$
$\frac{Z}{300}-\frac53x_1=x_2$
The slope is $-\frac53$
2. $\textrm{Second constraint}$
$6x_{2}\leq 240-10x_{1}$
$x_{2}\leq 40-\frac{10}6x_{1}$
The slope is $-\frac53$ as well.
So you may have multiple solution if at least one of the constraints have the same slope as the objective function. In this case you have two optimal BFS: $(x_1,x_2)=\left(\frac52, \frac{215}{6} \right), (x_1,x_2)=\left(15, 15 \right)$
And between the two corners on the green line the solution are optimal as well. The multiple optimal solution can be written as $(x_1^*,x_2^*)=\left(x_1, \ 40-\frac{10}6x_{1}\right) \ \ \forall \ \frac52 \leq x_1\leq 15$