We consider the problem of computing a bounded error approximation of the solution over a bounded time [0, T ], of a parameterized linear system, x(t) = Ax(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 eps > 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.