After receiving clarification on what you are trying to do through the chat linked hereto, here is the consolidated answer:

    myfit <- lm(X3 ~ X1 * X2, data=dt)
    newx <- model.matrix(myfit)
    solve(t(newx) %*% newx) %*% t(newx) %*% y

Essentially, you need to be sure that you have specified the correct parametrization for the model and correctly specified the model.  You can obtain the same parameterization as R using the "set first to zero" parameterization.  Take a look at the model matrix produced and stored in newx above.  This is how your model matrix must be specified in order to solve for $\hat{b}$ and replicate R's coefficients output.