Dedicated to Björn Engquist on the occasion of his 60th birthday.
Abstract.We explore an intermediate "Semi-Lagrangian" formulation of Geometric Optics type problems which stands between the classical "Lagrangian" Hamiltonian systems and its "Eulerian" Hamilton-Jacobi Partial Differential Equations counterpart. The goal is to design a numerical method which "trades" between advantages and drawbacks of both worlds to compute numerically the so-called "multi-valued" solution. We present a numerical algorithm and apply it to the "paraxial" simplification of 2-D geometric optics. We discuss its limitations and its possible extensions to higher dimensions. For an extended discussion on Lagrangian versus Eulerian methods in Geometric Optics and a review of some applications of these techniques see Engquist and Runborg, Acta Numer. (1 The Lagrangian/Eulerian setting. Our prototype problem is the following Hamiltonian system ẏ(s, y o ) = H p (s, y(s, y o ), p(s, y o )), y(0, y o ) = y o , p(s, y o ) = −H y (s, y(s, y o ), p(s, y o )), p(0, y o ) = p o (y o ), (1.1)where y o ∈ R d , s is a "time-like" dimension and( .) = d(.)ds . The notation f α will systematicaly denote partial differentiation of f with respect to α (and therefore a gradient when α ∈ R d with d > 1). The couples (y, p) are called "bicharacteristics" and lie in phase-space R s × R d y × R d p . The curve y(s, y o ) describes the trajectory of a ray issued from point y o and parameterized by s while p(s, y o )