I have read that when the regression model is in a linear form, using the OLS method is a good idea.
When instead, the model is not linear (for example, probit or logit) OLS method is not a good idea.
I am wondering, why? why does the OLS method work better when the model is linear?
As a general rule, fitting data that has an underlying/true model with that same model will usually produce the best results. More complex models may require more data to avoid overfitting, but that's not usually an issue with OLS vs. other methods when OLS describes the "true" model.
I think there should be clarification on what you are trying to ask, since the statements/assumptions in your OP aren't strictly true, although fitting a numeric response variable to OLS as a first model isn't necessarily a bad choice. If your response variable is just 1s and 0s, an OLS model is almost always a bad choice.
The OLS model is a way of fitting points to a line using a linear combination of X variables plus an estimate of normal error.
A linear model fits points to a line, although it might be transformed. In this case, a logistic regression with a logit/probit response function is considered a generalized linear model.
Alternatively, you could have a model that looks like
$y_i = \beta_0 + \beta_1 * x_1 + Laplace(0, b)$
where the error is described by the Laplacian distribution/double exponential function, or
$e^{-\left|{\frac{(y_i - \hat{y}_i)}{b}}\right|}$
In this case OLS would not be a good fit, despite it being a linear model, since the solution to this type of equation would require minimizing the sum of the absolute values of the errors, not the sum of the squares of the errors.
In most cases, though, you can get a decent idea of whether or not OLS is a good choice of models based on the shape of the distribution of the error residuals on a normal-normal plot.
e.g., compare the plots of normal error (left) with Laplacian error (right) below:
For a linear model the OLS estimator corresponds to the maximum-likelihood estimator (MLE), which has various good estimation properties. This is not true for non-linear models. In the latter case we can fit the model using the MLE or we can use iteratively reweighted least squares.
Minimizing square loss can be fine when the model is nonlinear, such as doing deep learning with mse as the loss function. Ordinary least squares, however, specifically refers to a linear model with the parameter vector estimated via $\hat{\beta} = (X^TX)^{-1}X^Ty$.
Thus, this is a bit of a tautology: we don't do ordinary least squares for nonlinear models because ordinary least squares is defined as applying to linear models.
It's totally fine to try out square loss in a logistic-style regression, however. In that situation, you would describe it as optimizing the so-called Brier score. This is not the standard way to do a logistic-style regression (which uses maximum likelihood, equivalent to minimizing crossentropy loss, not Brier score), but it might work well in some situations. I have even seen Yann LeCun post MNIST examples on his website, using neural netoworks, that optimize square loss (mse).
Be aware of how broad a linear regression can be, however. It has to do with linearity in the parameters, not just lines and higher-dimension analogues (e.g., planes). For instance, $y=\beta_0 + \beta_1x + \beta_2 x^2$ is a linear model, despite the parabolic curvature from the quadratic term.
