The saturation of I with respect J can be computed by computing until it stabilizes.
SINGULAR example:
ring R=...; ideal I=...; ideal J=...; int ii; I = std(I); while ( ii<=size(II)) { II=quotient(I,J); for ( ii=1; ii <=size(II); ii=ii++) { if (reduce(II[ii],I)!=0) break; } I=std(II); }// II is now (I:J).