At least the 2nd and 3rd solutions are correct. 

Your design matrix has dependent variables. For example the third column can be expressed in terms of the first two columns $x_3 = 0.15 + 0.45 x_2$ and the equation can also be expressed as 

$$\begin{array}{}
-0.2 + 0.4 x_ 1 + 1.2 x_3 &=& -0.2 + 0.4 x_ 1 + 1.2 (0.15 + 0.45 x_2) \\ &=& -0.2 + 0.4 x_ 1 + 1.2 (0.15 + 0.45 x_2)\\ &=& -0.02 + 0.94 x_1 \end{array}$$

### Methods 2 and 3

This last equation on the right hand side is the solution given by the 2nd and 3rd methods which probably drop one of the columns. In R you get the same behavior when we use the function `lm` which gives as output

    > lm(y~X+0)
    
    Call:
    lm(formula = y ~ X + 0)
    
    Coefficients:
       X1     X2     X3  
    -0.02   0.94     NA  

The last column is ignored when you give a computer the task to solve the equation.

### Method 1

Your 1st method probably attempts to inverse the (non-invertible) matrix anyway. For example, the inverse command does give some output. In my computer (an online https://www.tutorialspoint.com/execute_matlab_online.php) I get:

    disp(inv(A.'*A));
    -3.3777e+13  -1.0133e+14   2.2518e+14
    -1.0133e+14  -3.0399e+14   6.7554e+14
     2.2518e+14   6.7554e+14  -1.5012e+15
and gives some output that is close but not exact (possibly due to round of errors). 

In my case I got -0.071429  0.785714  0.342857, which is close to a correct solution $-0.071429+0.15 \cdot 0.342857 \approx -0.02$ and $0.785714+0.45 \cdot 0.342857 \approx 0.94$

In your case the difference is larger $-254.4 + 1696\cdot 0.15 \approx 0$ and $-762.3 + 1696*0.45 = 0.9$ (but this might be due to the output being given with less precision)

In R I can get the same result when I use the `solve` command while setting the tolerance parameter extremely low. In that case the inverse matrix is still computed; and it can be computed because the columns in the matrix X are not entirely dependent due to round off errors. 

    X = cbind(rep(1,20), seq(-1,1,length.out = 20), seq(-0.3,0.6,length.out = 20))
    beta = c(-0.2,0.4,1.2)
    y = X %*% beta
    X = round(X,5)
    solve(t(X) %*% X, tol = 10^-50) %*% t(X) %*% y
    #           [,1]
    #[1,] -0.0430674
    #[2,]  0.8707996
    #[3,]  0.1537806