This work studies the energy-critical inhomogeneous Schrödinger coupled equations with inverse square potential \(\begin{aligned} \text {i}\partial _t u_j +\Delta u_j-\frac{\lambda }{|x|^2}u_j =\pm |x|^{-\tau }\Big (\displaystyle \sum _{k=1}^ma_{jk}|u_k|^p\Big )|u_j|^{p-2}u_j. \end{aligned}\) Here and hereafter, \(u_j:\mathbb {R}\times \mathbb {R}^3\rightarrow \mathbb {C}\) and the above parameters satisfy \(1\le j\le m\) , \(\tau >0\) and \(\lambda >-\frac{1}{4}\) . The invariant Sobolev norm under the classical scaling \(\Vert \kappa ^\frac{2-\tau }{2(p-1)}{u_j}(\kappa ^{2} t,\kappa \cdot )\Vert _{\dot{H}^{s_c}}=\Vert {u_j}(\kappa ^{2} t)\Vert _{\dot{H}^{s_c}}\) gives the energy critical exponent \(p=3-\tau \) , which corresponds to \(1=s_c:=\frac{3}{2}-\frac{2-\tau }{2(p-1)}\) . In order to avoid a singular term \(|u_j|^{p-2}\) , one assumes that \(p\ge 2\) , which reads \(\tau <1\) . The assumption \(\lambda >-\frac{1}{4}\) is motivated by the critical Hardy inequality \(\frac{1}{4}\int _{\mathbb {R}^3}\frac{|f(x)|^2}{|x|^2}\,dx\le \int _{\mathbb {R}^3}|\nabla f(x)|^2\,dx\) , which guarantees that the operator \({\mathcal {K}}_\lambda :=-\Delta +\frac{\lambda }{|x|^2}\) is positive and gives the norm equivalence \(\Vert \nabla \cdot \Vert _{L^2(\mathbb {R}^3)}\simeq \Vert \sqrt{{\mathcal {K}}_\lambda }\cdot \Vert _{L^2(\mathbb {R}^3)}\) . The goal of this note is to develop a local theory and a global one in the energy space for small datum. One approaches with a classical fix point argument via two different Strichartz estimates. Indeed, in a first way, one uses some classical Schrödinger estimates via the fractional Hardy inequality \(\Vert |\cdot |^{-s}u\Vert _\beta \lesssim \Vert u\Vert _{\dot{W}^{s,\beta }}\) for \( 0<s<\frac{N}{\beta }\) and \(1<\beta <\infty \) , which enables to handle the inhomogeneous term \(|x|^{-\tau }\) in the source term. In a second method, one uses some adapted Strichartz estimates to weighted Sobolev spaces. This approach seems to be suitable to perform a finer analysis for the INLS model because the singularity in the nonlinear term can be handled more effectively in the weighted setting. This follows some ideas in Y. Lee and I. Seo (Arch. Math. (2021)). In both cases, one essential difference with the classical case \(\lambda =0\) is the need of some admissible pairs (q, r) satisfying the norm equivalence \(\Vert \nabla \cdot \Vert _{L^r(\mathbb {R}^3)}\simeq \Vert \sqrt{{\mathcal {K}}_\lambda }\cdot \Vert _{L^r(\mathbb {R}^3)}\) .