Hi Anomitra,
For the variable 'ne' you could define it as a Np x Np with diagonal elements as 0 and non diagonal elements as 1
ne= ne=ones(Np)-diag(ones(1,Np));
For H(x) you could you the function 'heaviside'
heaviside(x(j)+x(k)-xb(i))
Also to define B(i) you could use the for loop
B=zeros(1,Np)
for i=1:Np
for j=1:Np
for k=1:j
B(i)=B(i)+0.5*(1+ne(j,k))*heaviside(x(j)+x(k)+xb(i))-heaviside(x(j)+x(k)-xb(i))*(j+k)*N(j)*N(k);
end
end
end