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)
dt$X1
dt$X2
str(dt$X1)
attach(dt)
lm(X3~X1*X2)
newx<-model.matrix(lm(X3~X1*X2))
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.