How can I find a square root of $\mathbf A = \mathbf B^\top \mathbf B$, where $$ \mathbf B = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 \\ -r & 1 & 0 & \dots & 0 \\ 0 & -r & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & -r &1 \end{bmatrix} \in \mathrm{Mat}_{n\times n}(\mathbb R) $$
To be more precise, I have to estimate $\boldsymbol \beta$ from the following regression problem: $$ \mathcal N(\mathbf y \mid \mathbf X\boldsymbol \beta , \sigma ^2 \mathbf B^\top \mathbf B) $$ given $\mathbf y$ and $\mathbf X$. Let $\mathbf R = (\mathbf B^\top \mathbf B)^{\frac 12}$, then multiplying by $\mathbf R ^{-1}$: $$ \mathcal N(\mathbf R ^{-1}\mathbf y \mid \mathbf R ^{-1}\mathbf X\boldsymbol \beta , \sigma ^2 \mathbf I) $$ we obtain simple linear regression problem.
To find the square root of $A$ you can diagonalize it. In other words, you decompose $A = VDV^{-1}$, then compute $R= A^{1/2} = VD^{1/2}V^{-1}$. (Computing $D^{1/2}$ is fast since it's diagonal -- just the element-wise square root.)
Fortunately, $A$ is symmetric tridiagonal, so decomposing $A=VDV^{-1}$ is also computationally easy. (We have $A_{ii} = 1+r^2, A_{i+1,i}=A_{i,i+1}=-r \forall i\in[1,n-1]$ and $A_{nn}=1$.) If you store $A$ as a sparse matrix, Matlab should compute it quickly, and you could potentially try custom packages like trideigs for this specific case.
As an aside, for your problem you shouldn't actually invert $R$, but should just solve the system of equations.