I have a [data file][1] containing the four orbital element $(x,y,p_x,p_y)$ of a two-dimensional closed loop orbit. Note that the exact implicit equation of the orbit is not known. The plot of the 2D orbit on the $(x,y)$ plane data = Import["lyap_4d.dat", "Table"]; d00 = data[[All, {1, 2}]]; C0 = ListLinePlot[d00, AspectRatio -> 1, PlotRange -> All, Frame -> True, Axes -> False, PlotStyle -> {Black, Thick}] [![enter image description here][2]][2] Now let's define a rectangular grid of initial conditions around this orbit data = Flatten[Table[{i, j}, {i, 0.785, 1.135, 0.005}, {j, -0.45, 0.45, 0.005}], 1]; nic = Length[data] L0 = ListPlot[data, PlotStyle -> {Blue, PointSize[0.001]}]; plot1 = Show[{L0, C0}] [![enter image description here][3]][3] Finally, we determine which of the initial conditions are inside the orbit poly = Cases[Normal@C0, Line[x_] :> x, Infinity]; inPolyQ[poly_, pt_] := Graphics`Mesh`PointWindingNumber[poly, pt] =!= 0 data2 = Select[data, inPolyQ[poly[[1]], #] &]; nic2 = Length[data2] L1 = ListPlot[data2, PlotStyle -> {Red, PointSize[0.001]}]; plot2 = Show[{L1, C0}, Frame -> True, Axes -> False] [![enter image description here][4]][4] The list `data2` contains all the $(x_0,y_0)$ initial conditions inside the orbit. My question is the following: For the boundary orbit the $p_x$ and $p_y$ values are known. How can I use these values so as to predict in a way the $(p_{x0},p_{y0})$ of the orbits with initial conditions inside the orbit? In other words, how can I interpolate the $p_x$ and $p_y$ of the orbit by taking into account the values of the boundary? I use version 9.0 of Mathematica in Win XP SP3. Many thanks in advance! [1]: http://www.mediafire.com/download/cf1f0yz0jmdu1da/lyap_4d.dat [2]: https://i.sstatic.net/lMcQA.png [3]: https://i.sstatic.net/r5YHx.png [4]: https://i.sstatic.net/XhEWO.png