private solution after solving nonhomogenous euler equation 
Solve the equalation: $$x^2y''-3xy'+3y=\ln x$$

First I solved the homogeneous part and got $$y_h=c_1x+c_2x^3.$$ Then i wanted using variation of parameters writing that $$c_1' x'+c_2' (x^3)^\prime.$$ I solved the equations for $c_1,c_2$ and got a function which when i substitute in the equation, I get $x^2\ln x$ (not $\ln x$). I verified my answers with MuPaD and seems like the integration was correct. Am I using wrong method?

You have the correct homogeneous solution: $y_h = c_1 x+c_2x^3$.  
Now, let $y_1=x$ and $y_2=x^3$.  If we let $y_p=u_1(x)y_1(x) + u_2(x)y_2(x)$, our objective is to find $u_1(x)$ and $u_2(x)$ using Variation of Parameters.  When I took differential equations, we used Cramer's rule to show that if $y_p = u_1y_1 + u_2y_2$ was a solution to $y^{\prime\prime}+P(x)y^{\prime}+Q(x)y=f(x)$, then it followed that 
$$u_1^{\prime}(x) = \frac{\det\begin{bmatrix}0 & y_2\\ f(x) & y_2^{\prime}\end{bmatrix}}{\det\begin{bmatrix}y_1 & y_2\\ y_1^{\prime} & y_2^{\prime}\end{bmatrix}} = -\frac{y_2 f(x)}{W(y_1,y_2)}\tag{1}$$
and
$$u_2^{\prime}(x) = \frac{\det\begin{bmatrix}y_1 & 0\\ y_1^{\prime} & f(x)\end{bmatrix}}{\det\begin{bmatrix}y_1 & y_2\\ y_1^{\prime} & y_2^{\prime}\end{bmatrix}}= \frac{y_1f(x)}{W(y_1,y_2)}\tag{2}$$
where $W(y_1,y_2)=\det\begin{bmatrix}y_1 & y_2 \\ y_1^{\prime} & y_2^{\prime}\end{bmatrix}$ is the Wronskian of $y_1$ and $y_2$.  We also note that we can rewrite our original ODE as follows:
$$x^2y^{\prime\prime}-3xy^{\prime}+3y = \ln x \implies y^{\prime\prime} - \frac{3}{x}y^{\prime}+\frac{3}{x^2}y = \frac{\ln x}{x^2}.$$
Thus, $f(x)=\dfrac{\ln x}{x^2}$.  Plugging our $y_1$, $y_2$ and $f(x)$ into $(1)$ and $(2)$ gives us
$$u_1^{\prime}(x) = \frac{\det\begin{bmatrix}0 & x^3\\ \frac{\ln x}{x^2} & 3x^2\end{bmatrix}}{\det\begin{bmatrix}x & x^3\\ 1 & 3x^2\end{bmatrix}} = -\frac{\ln x}{2x^2}\qquad\text{ and } \qquad u_2^{\prime}(x) = \frac{\det\begin{bmatrix}x & 0\\ 1 & \frac{\ln x}{x^2}\end{bmatrix}}{\det\begin{bmatrix}x & x^3\\ 1 & 3x^2\end{bmatrix}} =\frac{\ln x}{2x^4}$$
We now integrate each one to see that $u_1(x)=\dfrac{1+\ln x}{2x}$  

 $$\begin{aligned}u_1(x) &= -\frac{1}{2}\int \frac{\ln x}{x^2}\,dx \\ &= -\frac{1}{2}\int te^{-t}\,dt;\quad(\text{sub: }t=\ln x)\\ &= -\frac{1}{2}\left[-te^{-t} + \int e^{-t}\,dt\right];\quad(\text{parts: }u=t,\,\,dv = e^{-t}\,dt)\\ &= \frac{1}{2}e^{-t}(1+t)\\ &= \frac{1+\ln x}{2x}\end{aligned}$$

and $u_2(x)=-\dfrac{1+3\ln x}{18x^3}$

 $$\begin{aligned}u_2(x) &= \frac{1}{2}\int\frac{\ln x}{x^4}\,dx \\ &= \frac{1}{2}\int te^{-3t}\,dt;\quad (\text{sub: } t=\ln x)\\ &= \frac{1}{2}\left[-\frac{1}{3}te^{-3t} + \frac{1}{3}\int e^{-3t}\,dt\right];\quad(\text{parts: }u=t,\,\,dv = e^{-3t}\,dt)\\ &= -\frac{1}{18}e^{-3t}(3t+1)\\ &= -\frac{1+3\ln x}{18x^3}\end{aligned} $$

Note that I didn't add any constants of integration to $u_1(x)$ and $u_2(x)$; the reason for this is because when we consider the general solution $y = y_h + y_p$, the constants of integration from this process get "absorbed" into the respective arbitrary coefficients of the homogeneous solution.
Therefore, the particular solution you should get is $$\begin{aligned}y_p &= u_1(x)y_1(x) + u_2(x)y_2(x) \\ &= \frac{1+\ln x}{2} - \frac{1+3\ln x}{18} \\ &= \frac{4}{9} +\frac{1}{3}\ln x \\ &= \frac{1}{9}(4+3\ln x)\end{aligned}$$
and hence the general solution is 
$$y=c_1x+c_2x^3+\frac{1}{9}(4+3\ln x)$$
which matches the solution given by WolframAlpha.

Alternatively, you could have avoided Variation of Parameters completely by solving this equation using the method of undetermined coefficients; in particular, you could have guessed a particular solution of the form $y_p = A+B\ln x$ and then you could have solved for $A$ and $B$.  You can read more about the method of undetermined coefficients for Cauchy-Euler equations here.
I hope this makes sense! 
