I have done this in the past for 3D hamiltonian systems in dynamical astronomy, so can broadly describe the adopted procedure.  Sorry if it sounds a bit technical.

1) The first step is to be able to find numerically the consequents of orbits in a Poincaré section using an ODE integrator.  Starting at position x in the section, one obtains the conseqent T(x).   Here it is important to get a precise consequent position.  Typically people tend to interpolate between successive time steps delivered by the ODE integrator when the section crossing has been detected, which is often not good enough.  An elegant method by Michel Hénon (*) is to integrate back to the section using the distance to the surface of section as independent variable, instead of time.  In one distance step the consequent T(x) is then found with numeric precision.   

2) The second step is to obain the linear map, a matrix, close to any consequent T(x).   Choosing an initial starting point x in the Poincaré section, one determines 2n consequents T(x+dz), T(x-dz), where dz are succesively all possible small variations from zero in only one of the n coordinates of the section.  The Jacobian matrix J is built from the n^2 numbers  (T(x+dz-T(x-dz)) with precision O(dz^2).

3) The third step is to use Newton method for finding the root of multidimensional functions.   What we want is to find the x* such that T(x*) = x*.  We get a better approximation of an initial guess x0 (Newton method) by

x1 = x0 - (J(x0) -I)^(-1).(T(x0)-x0)

It could be that the matrix J(x0) -I is badly conditionned. Instead of inverting J(x0) -I it is numerically better to solve the linear system 

(J(x0) -I)  Z = T(x0)-x0,

and then x1 = x0+Z

Once x1 is found one iterates to find x2, etc, until it converges to x* to the wished level of precision.   The convergence is quadratic, at each iterations the accurate digits double. 

4) This method yields the stable or unstable periodic solutions of T.  The degree of stability or instability near x* is given by the eigenvalues of the final J(x*).   Any good eigenvalue finder can yields then real or complex eigenvalues.  If any real part of the eigenvalues are positive, the periodic orbit x* is unstable.

 (*) Hénon, M.:  Physica D: Nonlinear Phenomena, Volume 5, Issue 2-3, p. 412-414.  (1982)

Go in reverse time in the numerical method for attraction to the periodic orbit (if unstable periodic orbit is not the saddle type). You see the articles "A new accurate numerical method of approximation of chaotic solutions of dynamical model equations with quadratic nonlinearities" and "A new reliable numerical method for computing chaotic solutions of dynamical systems: the Chen attractor case".

Article A New Reliable Numerical Method for Computing Chaotic Soluti...

Article A New Accurate Numerical Method of Approximation of Chaotic ...

One simple solution to try is to integrate backwards. Not always works, I guess.

A more systematic approach, I think, would be to extend your ode to an evolution PDE.

For example, suppose you have a driven-damped oscillator with quadratic nonlinearity  written for its coordinate f,  and where x-variable, is your 'time'. 

Just add a time derivative $\partial f/\partial t$ with a proper sign and you'll  get a KdV-Burgers equation which you can always solve numerically, or even analytically in certain cases. If e.g., the driver is periodic, you'll easily get a periodic solution,  period-doubling bifurcation, and even chaos in x. All those will be stable attractors of the extended evolution equation (KdVB). See a paper attached to this.

Article Spatial chaos in weakly dispersive and viscous media: A nonp...

There is a simple way to extract your periodic orbits from your previous work. When you write "I could get the periodic points through first return map", I assume you obtained the points in the return maps. For instance for a period-2 orbits, you have the coordinate in the format (x_n , x_n+1) for the two points of the orbits:

(0.3, 0.5)

(0.5, 0.3)

When you extract these points, you can also collect the iteration step of your solver:

(0.3, 0.5) --- step k

(0.5, 0.3) --- step k+j

As your system is deterministic, you can run the same simulation solver and collect the j points (x, y, z) of your period-2 orbits between the steps k and k+j.

It works for all types of orbits. I use it for several systems such as Rössler, Lorenz and Chen systems.

First of all, backward integration does not help since the orbits in a chaotic attractor are unstable forward and backward. (A heuristic argument: your system has at least one  positive Lyapunov exponent since it is chaotic, and it has at least one negative exponent since it is an attractor where the sum of all exponents is negative. Reversing time reverses sign of the exponents and you get a system having now negative and positive exponents.)

Any software system dedicated to numerical analysis of dynamical systems is able to compute periodic orbits independent of their stability (but stability affects how many of the infinitely many periodic orbits can be obtained in practice). These software system do professionally what Daniel Pfenniger proposed to do.

When I was working in this field (about 20 years ago) corresponding free software systems were AUTO, CONTENT (I think both are still alive), and my own CANDYS/QA (no longer supported). Maybe that nowadays also MATLAB supports computing of unstable orbits.

Wolfgang, if the limiting trajectory has the saddle type, then it is unstable in both directions in time. The backward integration with high accuracy in this case can give a check of results obtained forward in time.

Alexander, you are right, backward integration may serve as check for forward integration. But you don't get any periodic orbits this way (and plotting embedded periodic orbits were the original problem). .

Wolfgang, at the beginning we move the given initial point by the corresponding trajectory on a long time for better approximation to the attractor. Further we research the return of trajectory in a neighborhood of the obtained initial point because the limiting trajectory on the attractor has a property of recurrence, i.e. we need to explore a rapprochement of the trajectory with the initial point and to fix the moments of time when it occurs.

Alexander, approaching an orbit's starting point by recurrence is not yet a periodic orbit, it is merely a possible initial point (and recurrence time is an approximate period length, both are needed) for a Newton-like method to actually compute a periodic orbit. This is what the mentioned software systems do: computing an periodic orbit independent of its stability from a reasonable approximation. The difference between start and end points of such an orbit is not zero because integration errors etc. but it is much smaller than the difference you generally get by recurrence. Sure, you may feed the point and time you found by recurrence as initial approximation into the software system (depending on the software, such an approximation may be requested). For example, CANDYS/QA had the option to look for a reasonable recurrence point as initial approximation (what you propose was there a preparation step for the proper work).

I agree with you but at least we can be judged by the distribution of times of the return points what is the approximate behavior of recurrent trajectory: if the moments of time of return of trajectory to the neighborhood of the initial point are approximately at the same distance from each other, then the corresponding trajectory describes by the almost periodic function (the periodic solution is a particular case). Of course, this argument is some approximation (we can not calculate indefinitely).

I've had good success with a simple, but not very efficient method. Using a small but fixed time step size, I save the time series of the state variables of the N most recent points on the orbit in a circular buffer. When each new point is added to the buffer, I calculate its Euclidean distance from each point in the buffer and save the value of the smallest separation. When the orbit eventually comes sufficiently close to a UPO as indicated by an acceptably small separation, I then just plot the orbit from the points stored in the buffer. The position of the points in the buffer gives the period of the orbit to the precision of the time step size. Since a chaotic system will have infinitely many UPOs mostly with very long periods, you should use a value of N small enough to capture only the ones with small periods that presumably most interest you. I usually start with a small N and slowly increase its value to identify a succession of UPOs with progressively longer periods. I have a few other tricks, but this basic method usually suffices if you are patient and do not need high precision. A similar method is described in more detail along with examples from the Lorenz attractor in my "Chaos and Time-Series Analysis" book (Oxford, 2003).

You could look at Cell Mapping methods, like the Interpolated Cell Mapping. A starting point Could be Tongue & Gu (1988)

Article A Theoretical Basis for Interpolated Cell Mapping