Wrong descent direction in ISTA Algorithm

The issue is about this paper cited as : arXiv:1309.6683 (Dynamic Structural Equation Models for Social Network Topology Inference)

In this paper I want to compute the ISTA Algorithm using MATLAB (Pseudo real-time ISTA for topology tracking), I followed all the steps as cited in the paper, but I am getting wrong results. I tested whether the descent direction is correct or not, and it was wrong (according to my test in which I plot the cost function w.r.t. the step size).

 % Test the descent direction (wrt A)

                stepsize_vec=linspace(0,1,1000);

                fct=zeros(1,1000);

                for ip=1:1000

                    % Update_of_A

% Soft_thresholding function gets as inputs : row a_i of A , the Lipschitz constant and the gradient_a_i 

                    Update_of_A(ii,ind_vec_ii)=Soft_thresholding(transpose(A(ii,ind_vec_ii)),(1/ip),gradient_a_i);

                    fct(ip)=norm(Y_t - Update_of_A*Y_t - B*Y_t ,'fro')^2;

                end

                figure;

                xlabel('stepsize')

                ylabel('cost function')

                plot(stepsize_vec,fct)


      % Test the descent direction (wrt B)

%                 stepsize_vec=linspace(0,1,1000);

%                 fct=zeros(1,1000);

%                 for ip=1:1000

%                     % Update_of_B

%                     Update_of_B(ii,ind_vec_ii)=B(ii,ii) - (gradient_bii/(1/ip));

%                     fct(ip)=norm(Y_t - A*Y_t - Update_of_B*Y_t ,'fro')^2 ;

%                 end

%                 figure;

%                 xlabel('stepsize')

%                 ylabel('cost function')

%                 plot(stepsize_vec,fct)

I am getting monotonically decreasing functions.

The gradients are computed as follows :

gradient_a_i = Sigma_t(ind_vec_ii,ind_vec_ii) * A(ii,ind_vec_ii).' + Ybarre_t(ind_vec_ii,:) * X(ii,:).' * B(ii,ii) - Sigma_t(ind_vec_ii,ii);            

gradient_bii=A(ii,ind_vec_ii) * Ybarre_t(ind_vec_ii,:) * X(ii,:).' + ( ((1-beta.^t)/(1-beta))*B(ii,ii)* X(ii,:)*X(ii,:).’ )- Ybarre_t(ii,:) * X(ii,:).’;

I don't know if I'm getting wrong results because my test is not correct or the function it self is wrong. Do you think that the gradient is computed in a wrong way ?? Or what could be the reason why I am getting wrong results (wrong descent direction) ?