>Is it correct to say that the non-linear activation function's main purpose is to allow the neural network's decision boundary to be non-linear?

Yes.

Neural networks compose several functions in layers: the output of a previous layer is the input to the next layer. If you compose linear functions, these functions are all linear. So the result of stacking several linear functions together is a linear function.  Showing this is simple algebra:

$$
\begin{align}
\hat{y} &= W_2(W_1x + b_1)+b_2 \\
&=  \underbrace{W_2W_1}_W x+\underbrace{W_2b_1+b_2}_b \\
&= Wx+b
\end{align}
$$


Any model which minimizes a loss $L(y,\hat{y})$ over parameters $W_1,W_2,b_1,b_2$ is equivalent to a model which minimizes the same loss over parameters $W,b$. In the case that the loss is the square error loss, this is exactly the same as an OLS model. 

On the other hand, using a nonlinear function makes the map from the input to the output nonlinear.$$
\hat{y} = f(W_2 f(W_1x + b_1)+b_2) \\
$$
 For some activation function $f$, such as $\tanh$ or ReLU, this cannot be rewritten as a single linear operation on $x$.

If you are estimating this model and observe a discrepancy between an OLS solution and a neural network optimized using gradient descent, it's probably due to either or both of two facts (1) gradient descent is not an effective optimizer for certain problems; (2) the problem is ill-conditioned. For more information, see 

- [Reflections on Random Kitchen Sinks][1]
- https://stats.stackexchange.com/questions/160179/do-we-need-gradient-descent-to-find-the-coefficients-of-a-linear-regression-mode/164164#164164

- https://stats.stackexchange.com/questions/278755/why-use-gradient-descent-for-linear-regression-when-a-closed-form-math-solution/278794#278794

This isn't unique to classification problems, either. If you have some sort of regression problem (such as an output that can take on any real number), then using nonlinear activation functions is necessary to model a nonlinear relationship between the inputs and outputs.

For example, a ReLU function's output is either 0 or positive. If the unit is 0, it is effectively "off," so the inputs to the unit are not propagated forward from that function. If the unit is on, the input data is reflected in subsequent layers through that unit. ReLU itself is not linear, and neither is the composition of several layers of several ReLU functions. So the mapping from inputs to classification outcomes is not linear either.


  [1]: http://www.argmin.net/2017/12/05/kitchen-sinks/