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.

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.

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.

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.