For solving this problem you can use the method of biconjugate gradient.

For example see: http://astronu.jinr.ru/wiki/upload/d/d6/NumericalRecipesinC.pdf

I have used this method but I got the wrong answer .

I also have used this method and I have obtained good results.

The book "NumericalRecipesinC" is very useful because it include C++ codes

Also you can see next books:

1)Bunch, J.R., and Rose, D.J. (eds.) Sparse Matrix Computations. New York: AcademicPress.1976

2)Duff, I.S., and Stewart, G.W. (eds.) 1979,Sparse Matrix Proceedings Philadelphia:S.I.A.M., 1978

use sundials framework

I do emphesize, I need a subroutine, not a software, package or library !

Sundials library has many subroutins that you can customize to solve your problem

How should I implement it ?

Is there any free web link ?

Even if you are using C++, your sparse matrix will consist of three arrays (one for the data, one for the column indices, and one for the size of the rows), which you could extract and feed even to C routines/libraries. If it uses std::vector for the data just use &vector[0] to convert to C arrays. Thus you could use SuperLU, MUMPS or others to do the actual solve

you can try LU method in the Direct Methods for Sparse Linear Systems

https://www.cise.ufl.edu/research/sparse/CSparse