We study the existence of analytic solutions of the system of difference equations given by y(epsilon, x + epsilon) = F(epsilon, x,a(epsilon), y(epsilon, x)), where F is an analytic function of epsilon, x, a and y on a certain neighborhood of a point (0, x(0), a(0), y(0)) is an element of C x C x C-m x C-n with y(0) = F(0, x(0), a(0), y(0)). We assume the existence of a "slow curve" phi(0)(x) satisfying phi(0)(x) = F(0, x, a0, phi(0)(x)) in a neighborhood of x(0) and phi(0)(x(0)) = y(0). Under the assumption that the Jacobian (partial derivative F/partial derivative y)(0, x, a(0), phi(0)(x)) - I is invertible except at x = x(0), we first show, under some transversality condition, the existence of a unique formal solution ((a) over cap(epsilon)(y) over cap(epsilon, x)) and establish its Gevrey-1 character. Then, we construct a quasi-solution, ((a) over tilde(epsilon), (y) over tilde(epsilon, x)) in the sense that it solves our system equation up to an exponentially small remainder. Next, we show the existence of an actual analytic solution (a(epsilon), y(epsilon, x)), defined for epsilon in some sector with vertex at epsilon = 0 and in some neighborhood of x = x(0) ( to be constructed), that is exponentially close to the above quasi-solution. We also show the exponential closeness of any two quasi-solutions ( and hence also of actual solutions) tending to the same slow curve. Finally, as an application of this theory, we give two numerical examples with some figures of discrete bounded solutions for two difference equations of second order.