Solving inhomogenous ODE I have an inhomogenous  ODE. The main issue here is variables are matrices. It is bit of matrix calculus.    A solution would be highly appreciated interms of x . I guess we can use same methods for solving ODEs but have to be careful because these are matrices
$R'(x)-(C_1 +C_2 x) R(x) = R_1-C_1 R_0\, x $ 
where  except x rest are $3\times 3$ matrices means $C_1,C_2,R'(x),R(x)$  all are matrices. x is a scalar variable. $ C_1,C_2,R_0 $ are constant $3\times 3$  matrices . .$C_1$ and $C_2$ are skew symmetric matrices   
 A: My answer for this will overlap quite a bit with my answer for your earlier question but I'll try to elaborate on the more specifically ODE and matrix aspects. My approach will look a bit different than @RRL's, but the content is essentially the same (and I'm choosing my notation to emphasize that).
First, some notation and a review of the scalar case. Note that may write $$R'(x)-A(x) R(x) = G(x)$$ where we have defined $A(x) = C_1+C_2 x$ and $G(x) = R_1 - C_0 x$. In the scalar case, we could then observe straightforwardly observe that the LHS has an integrating factor $\Phi(x)=\exp\left[-\int_0^{x} A(x') dx'\right]$. (In your case, $\Phi(x)=e^{-C_1 x-C_2 x^2/2}$; I have taken the lower bound to match the boundary condition at $x=0$ from the original problem). What that means is we may write 
$$ \frac{d}{dx} \Phi(x) R(x) = \Phi(x)R'(x)+\Phi'(x) R(x) = \Phi(x)\left(R'(x)-A(x)R(x)\right)=\Phi(x) G(x)\hspace{.4cm}\text{(1)}$$ where we have exploited the fact that $\Phi'(x) = -\Phi(x) A(x)$ from the fundamental theorem of calculus. We can then integrate both sides from $0$ to $x$ and rearrange both sides to obtain the formal solution
$$R(x)=\Phi^{-1}(x)\left[R(0) + \int_0^{x} \Phi(x')G(x')\,dx'\right]  =\Phi^{-1}(x)R(0) + \int_0^{x} \Phi(x'-x)G(x')\,dx'\hspace{.4cm}\text{(2)}$$ where in the second equality we have taken advantage of $\Phi(x)$ being the exponent of a definite integral. So we've converted  the problem of finding $R(x)$ to computing a particular integral.
The question is, what changes when we go to matrices in the coefficients? An obvious issue is that $\Phi(x)$ now will be the exponential of a matrix-valued function. This isn't a knock-down problem, though: define the matrix exponential as the infinite sum $\exp(A) =\sum_{k=0}^\infty A^k/k!$, and observe that one still has the property $$\frac{d}{dx} \exp{(x A)}=\frac{d}{dx}\sum_{k=0}^\infty A^k \frac{x^k}{k!}=\sum_{k=1}^\infty A^k \frac{x^{k-1}}{(k-1)!} = A\cdot \exp(x A).$$ As a consequence, the $\Phi(x)$ defined earlier has the same derivative and so equation (1) is still valid. The same moreover holds for equation (2), with $\Phi(x)^{-1}$ understood to be the inverse matrix (the second equality is now valid only if $A(x)$ commutes with $A(x')$).
