Ans: Let’s start by defining matrix Q.

library(matlib)
Q <- matrix(c(0, 2/3, 0, 1/3, 0, 2/3, 0, 1/3, 0),nrow = 3, byrow = T)
Q
##           [,1]      [,2]      [,3]
## [1,] 0.0000000 0.6666667 0.0000000
## [2,] 0.3333333 0.0000000 0.6666667
## [3,] 0.0000000 0.3333333 0.0000000
I <- diag(1,3,3)
I
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

Calculating for N,

N <- inv(I - Q)
N
##      [,1] [,2] [,3]
## [1,]  1.4  1.2  0.8
## [2,]  0.6  1.8  1.2
## [3,]  0.2  0.6  1.4

Calculating for Nc,

c <- matrix(1, nrow = 3)
Nc <- N %*% c
Nc
##      [,1]
## [1,]  3.4
## [2,]  3.6
## [3,]  2.2

Calculating for NR,

# By looking at the definition in the canonical form, we can deduce that R is as follows
R <- matrix(c(1/3,0,0,0,0,2/3), byrow = T, nrow = 3)
R
##           [,1]      [,2]
## [1,] 0.3333333 0.0000000
## [2,] 0.0000000 0.0000000
## [3,] 0.0000000 0.6666667
NR <- N %*% R
NR
##            [,1]      [,2]
## [1,] 0.46666667 0.5333333
## [2,] 0.20000000 0.8000000
## [3,] 0.06666667 0.9333333

So my solution set is:

\(N = \begin{pmatrix}1.4 & 1.2 & 0.8\\ 0.6 & 1.8 & 1.2 \\ 0.2 & 0.6 & 1.4 \end{pmatrix}\)

Nc\(=\begin{pmatrix}3.4\\ 3.6\\ 2.2\end{pmatrix}\)

NR = \(\begin{pmatrix}4.67 & 0.53\\ 0.20 & 0.80\\ 0.07 & 0.93\end{pmatrix}\)