We consider the problem of computing a bounded error ap-proximation of the solution over a bounded time [0, T], of a parameterized linear system, ẋ(t) = A x(t), where A is constrained by a compact polyhedron Ω. Our method consists of sampling the time domain [0, T] as well as the parameter space Ω and constructing a continuous piecewise bilinear function which interpolates the solution of the parameterized system at these sample points. More precisely, given an ϵ > 0 , we compute a sampling interval δ > 0, such that the piecewise bilinear function obtained from the sample points is within of the original trajectory. We present experimental results which suggest that our method is scalable.