diff --git a/pub/lib/olf.cpp b/pub/lib/olf.cpp
index 7e01b2d9427d4f644abeb12f3210cf1f3385f660..80138272f1e13849c51dc569f375345f4c869717 100644
--- a/pub/lib/olf.cpp
+++ b/pub/lib/olf.cpp
@@ -387,6 +387,19 @@ int COld::nPts() const
     return np;
 }
 
+//! Return maximum number of points.
+
+int COld::nPtsMax() const
+{
+    if (!nJ())
+        throw "BUG: nPtsMax() called while nSpec=0";
+    int np = VS(0)->size();
+    for (int j = 0; j < nJ(); ++j)
+        if (VS(j)->size() > np)
+            np = VS(j)->size();
+    return np;
+}
+
 
 //! Returns true if all data points have nonzero dy.
 
diff --git a/pub/lib/olf.hpp b/pub/lib/olf.hpp
index 55c3a38cf2662826d849bee9f515adc1ec8d2761..1ffee0896fa64bd40181c3d8fb95fb73f44c3013 100644
--- a/pub/lib/olf.hpp
+++ b/pub/lib/olf.hpp
@@ -79,6 +79,7 @@ class COld : public COlo {
 
     int nPts(int j) const;
     int nPts() const; // 0 unless all spectra j have same nPts(j)
+    int nPtsMax() const;
     bool dy_are_all_nonzero() const;
     bool has_dy() const;
     bool x_sorted() const;
diff --git a/pub/lib/opr.cpp b/pub/lib/opr.cpp
index 1a719b40e377a511b667428fd0f70d70c49bc52a..3c813fbb75371ddc5cd821c45ae2551b33b653bd 100644
--- a/pub/lib/opr.cpp
+++ b/pub/lib/opr.cpp
@@ -345,61 +345,86 @@ void NOperate::IntXY(string mode)
     else
         throw S("invalid command");
 
-    static int icolx = 0, icoly = 1, icold = 2;
+    static int icolx = 0;
+    static string sely, seld;
     static CCoord xco, yco;
+
     icolx = iask("x from column", icolx);
     if (icolx < 0)
         return;
-    icoly = iask("y from column", icoly);
-    if (icoly < 0)
-        return;
+
+    sely = wask("y from column(s)", sely);
+    if (sely == "")
+	    return;
+    RNode Tsely = user_xaxparse(sely.c_str());
+
+    RNode Tseld;
     if (with_d) {
-        icold = iask("dy from column", icold);
-        if (icold < 0)
-            return;
+	seld = wask("dy from column(s)", seld);
+	if (seld == "")
+	    return;
+	Tseld = user_xaxparse(seld.c_str());
     }
 
     while (const COld* fin = fiter.nextD()) {
         int nz = fin->ZCo.size();
 
-        POld fout(fin->new_POld());
-        string lin = S("oi") + mode + " " + S(icolx) + " " + S(icoly);
-        if (with_d)
-            lin += " " + S(icold);
-        fout->log_action(lin);
-
-        fout->xco = CCoord("x", "");
-        fout->yco = CCoord("y", "");
-        if (nz)
-            fout->ZCo.pop_back();
-
-        PSpec sout(new CSpec);
-        sout->z = fin->V[0]->z;
-        if (nz)
-            sout->z.pop_back();
-
-        for (int j = 0; j < fin->nJ(); j++) {
-            int n = fin->nPts(j);
-            if (icolx >= n || icoly >= n || (with_d && icold >= n))
-                throw "not enough columns in spectrum " + S(j);
-            if (with_d)
-                sout->push_xyd(fin->VS(j)->y[icolx], fin->VS(j)->y[icoly], fin->VS(j)->y[icold]);
-            else
-                sout->push_xy(fin->VS(j)->y[icolx], fin->VS(j)->y[icoly]);
-            sout->z = fin->VS(j)->z;
-            if (nz > 1) { // new spectrum if jump in other z value
-                if (j + 1 < fin->nJ()
-                    && fin->V[j + 1]->z[nz - 2]->to_r() != fin->V[j]->z[nz - 2]->to_r()) {
-                    fout->V.push_back(move(sout));
-                    sout = PSpec(new CSpec);
-                    sout->z = fin->V[j + 1]->z;
-                    if (nz)
-                        sout->z.pop_back();
-                }
-            }
-        }
-        fout->V.push_back(move(sout));
-        SMem::instance()->mem_store(move(fout), fiter.k());
+	RObjVecInt vsely = Tsely->to_index_list(fin->nPtsMax());
+	if (!vsely)
+	    throw S("not a valid index list for y");
+	int nfout = vsely->size();
+	RObjVecInt vseld;
+	if (with_d) {
+	    vseld = Tseld->to_index_list(fin->nPtsMax());
+	    if (!vseld)
+		throw S("not a valid index list for dy");
+	    if (vseld->size() != nfout)
+		throw S("index lists for y and dy have different sizes");
+	}
+
+	for (int ifout = 0; ifout < nfout; ++ifout) {
+	    const int icoly = vsely->to_i(ifout);
+	    const int icold = with_d ? vseld->to_i(ifout) : -1;
+
+	    POld fout(fin->new_POld());
+	    string lin = S("oi") + mode + " " + S(icolx) + " " + S(icoly);
+	    if (with_d)
+		lin += " " + S(icold);
+	    fout->log_action(lin);
+
+	    fout->xco = CCoord("col"+S(icolx), "");
+	    fout->yco = CCoord("col"+S(icoly), "");
+	    if (nz)
+		fout->ZCo.pop_back();
+
+	    PSpec sout(new CSpec);
+	    sout->z = fin->V[0]->z;
+	    if (nz)
+		sout->z.pop_back();
+
+	    for (int j = 0; j < fin->nJ(); j++) {
+		int n = fin->nPts(j);
+		if (icolx >= n || icoly >= n || (with_d && icold >= n))
+		    throw "not enough columns in spectrum " + S(j);
+		if (with_d)
+		    sout->push_xyd(fin->VS(j)->y[icolx], fin->VS(j)->y[icoly], fin->VS(j)->y[icold]);
+		else
+		    sout->push_xy(fin->VS(j)->y[icolx], fin->VS(j)->y[icoly]);
+		sout->z = fin->VS(j)->z;
+		if (nz > 1) { // new spectrum if jump in other z value
+		    if (j + 1 < fin->nJ()
+			&& fin->V[j + 1]->z[nz - 2]->to_r() != fin->V[j]->z[nz - 2]->to_r()) {
+			fout->V.push_back(move(sout));
+			sout = PSpec(new CSpec);
+			sout->z = fin->V[j + 1]->z;
+			if (nz)
+			    sout->z.pop_back();
+		    }
+		}
+	    }
+	    fout->V.push_back(move(sout));
+	    SMem::instance()->mem_store(move(fout), fiter.k());
+	}
     }
 }