A general procedure, based on the knowledge of the input and the output signals, is proposed to approximate the prescribed linear time-invariant (LTI) systems by means of optimal Laguerre models. The main contribution of this paper is to apply the Newton Raphson＇s iterative technique to compute the so-called optimal Laguerre pole in a continuous-time case (or optimal time scale factor in a discrete-time case) and especially to show that the gradient and the Hessian can be expressed analytically. Moreover, the excitations used are not limited to the ones that ensure the orthogonality of the outputs of Laguerre filters (i.e., Dirac delta or white noise) as is usually done in existing methods, however persislently exciting input signal(s) are used. The proposed procedure will be directly formulated for multi-input/single-output (MISO) systems, single-input/single-output (SISO) systems being a special case with the number of inputs equal to one. The proposed algorithm has direct applications in system identification, model reduction, and noisy modeling.