Overview
I don't think that solving classification problems using linear regression is usually the best approach (see notes below), but it can be done. For multiclass problems, multinomial logistic regression would typically be used rather than a combination of multiple regular logistic regression models. By analogy, one could instead use least squares linear regression with multiple outputs.
Approach
Suppose we have training data $\big\{ (x_i,y_i) \big\}_{i=1}^n$ where each $x_i \in \mathbb{R}^d$ is an input point with class label $y_i$. Say there are $k$ classes. We can represent each label as a binary vector $y_i \in \{0,1\}^{k}$, whose $j$th entry is $1$ if point $i$ is a member of class $j$, otherwise $0$. The regression problem is to predict the vector-valued class labels as a linear function of the inputs, such that the squared error is minimized:
$$\min_W \ \sum_{i=1}^n \|y_i - W x_i\|^2$$
where $W \in \mathbb{R}^{k \times d}$ is a weight matrix and $\|\cdot\|^2$ is the squared $\ell_2$ norm. The inputs should contain a constant feature (i.e. one element of $x_i$ should always be $1$), so we don't have to worry about extra bias/intercept terms.
To predict the class for a new input $x$, compute the vector $a = W x$, where $a_i$ is the projection of the input onto the the $i$th row of $W$ (the weights for the $i$th class). Then, some rule can be applied to map the projections to a single class. For example, we could choose the class with the maximal projection: $\arg \max_i a_i$. This is loosely analogous to selecting the most probable class in multinomial logistic regression.
Example
Here's a plot of the decision boundaries learned from a set of 2d points, using the above method. Colors represent true class labels.
Notes
This method sacrifices the principled, probabilistic approach used in multinomial logistic regression. The squared error is also an odd choice for classification problems, where we're predicting binary values (or binary vectors, as above). The issue is that the squared error penalizes large outputs, even when these ought to be considered correct. For example, suppose the true class label is $[1,0,0]$. Outputting $[2,0,0]$ (which should correspond to high confidence in the correct class) is just as costly as outputting $[0,0,1]$ (which corresponds to high confidence in the wrong class). Even if one is willing to abandon probabilistic models, there are other loss functions designed specifically for classification, like the hinge loss used in support vector machines. The main benefit of the squared error is computational efficiency. But, this doesn't seem particularly necessary in most cases, given that we can routinely solve much more complicated problems involving massive datasets. Nevertheless, one does sometimes see the squared error used in the literature for classification problems (apparently with success). Least squares support vector machines are the most prominent example that comes to mind.
Code
Matlab code to generate the example plot above. Matrices are transposed relative to the text above, since points and labels are stored as rows.
%% generate toy dataset
% how many points and classes
n = 300;
k = 3;
% randomly choose class labels (integers from 1 to k)
c = randi(k, n, 1);
% convert labels to binary indicator vectors
% Y(i,j) = 1 if point i in class j, else 0
Y = full(sparse((1:n)', c, 1));
% mean of input points in each class
mu = [
0, 0;
4, 0;
0, 4
];
% sample 2d input points from gaussian distributions
% w/ class-specific means
X = randn(n, 2) + mu(c, :);
% add a column of ones
X = [X, ones(n,1)];
%% fit weights using least squares
W = X \ Y;
%% out-of-sample prediction
% generate new test points on a grid covering the training points
[xtest2, xtest1] = ndgrid( ...
linspace(min(X(:,2)), max(X(:,2)), 501), ...
linspace(min(X(:,1)), max(X(:,1)), 501) ...
);
X_test = [xtest1(:), xtest2(:)];
% add a column of ones
X_test = [X_test, ones(size(X_test,1), 1)];
% project test points onto weights
A_test = X_test * W;
% predict class for each test point
% choose class w/ maximal projection
[~, c_test] = max(A_test, [], 2);
%% plot
% plot decision boundary
% using contour plot of predicted class labels at grid points
figure;
contour(xtest1, xtest2, reshape(c_test, size(xtest1)), 'color', 'k');
% plot training data colored by true class label
hold on;
scatter(X(:,1), X(:,2), [], c, 'filled');