  In this paper, we investigative the large time decay and stability to any
given global smooth solutions of the $3$D incompressible inhomogeneous MHD
systems. We proved that given a solution $(a, u, B)$ of (\ref{mhd_a}), the
velocity field and magnetic field decay to $0$ with an explicit rate, for $u$
which coincide with incompressible inhomogeneous Navier-Stokes equations
\cite{zhangping}. In particular, we give the decay rate of higher order
derivatives of $u$ and $B$ which is useful to prove our main stability result.
For a large solutions of (\ref{mhd_a}) denoted by $(a, u, B)$, we proved that a
small perturbation to the initial data still generates a unique global smooth
solution and the smooth solution keeps close to the reference solution $(a, u,
B)$. Due to the coupling between $u$ and $B$, we used elliptic estimates to get
$\|(u, B)\|_{L^{1}(\mathbb{R}^{+};\dot{B}_{2,1}^{5/2})} < C$, which is
different to Navier-Stokes equations.
