The Szekeres inhomogeneous models can be used to model the true lumpy universe that we observe. This family of exact solutions to Einstein's equations was originally derived with a general metric that has no symmetries. In this work, we develop and use a framework to integrate the angular diameter and luminosity distances in the general Szekeres models. We use the affine null geodesic equations in order to derive a set of first-order ordinary differential equations that can be integrated numerically to calculate the partial derivatives of the null vector components. These equations allow the integration in all generality of the distances in the Szekeres models and some examples are given. The redshift is determined from simultaneous integration of the null geodesic equations. This work does not assume spherical or axial symmetry, and the results will be useful for comparisons of the general Szekeres inhomogeneous models to current and future cosmological data.

We use the Szekeres inhomogeneous relativistic models in order to fit supernova combined data sets. We show that with a choice of the spatial curvature function that is guided by current observations, the models fit the supernova data almost as well as the LCDM model without requiring a dark energy component. The Szekeres models were originally derived as an exact solution to Einstein's equations with a general metric that has no symmetries and are regarded as good candidates to model the true lumpy universe that we observe. The null geodesics in these models are not radial. The best fit model found is also consistent with the requirement of spatial flatness at CMB scales. The first results presented here seem to encourage further investigations of apparent acceleration using various inhomogeneous models and other constraints from CMB and large structure need to be explored next.