function [up,rhs] = boundary(neq,jdiag,value,ic,up,rhs)

for j=1:neq
    if (jdiag(j,1) < jdiag(ic,1))
        k1 = ic - j;
        k2 = jdiag(ic,1) - jdiag(ic-1,1) - 1;
        if (k1 <= k2)
            rhs(j,1) = rhs(j,1) - up(jdiag(ic,1)-k1,1)*value;
	        up(jdiag(ic)-k1,1) = 0;
        end
    elseif (jdiag(j,1) == jdiag(ic,1))
        rhs(ic,1) = value;
        up(jdiag(j,1),1) = 1;
    else
        k1 = j - ic;
	    k2 = jdiag(j,1) - jdiag(j-1,1) - 1;
	    if (k1 <= k2)
            rhs(j,1) = rhs(j,1) - up(jdiag(j,1)-k1,1)*value;
	        up(jdiag(j)-k1,1) = 0;
        end
    end
end