The numerical parametrization method (PM), originally created for optimal control problems, is specificated for classical calculus of variation problems that arise in connection with singular implicit (IDEs) and differential-algebraic equations (DAEs). The PM for IDEs is based on representation of the required solution as a spline with moving knots and on minimization of the discrepancy functional with respect to the spline parameters. Such splines are named variational splines. For DAEs only finite entering functions can be represented by splines, and the functional under minimization is the discrepancy of the algebraic subsystem. The first and the second derivatives of the functionals are calculated in two ways -for DAEs with the help of adjoint variables, and for IDE directly. The PM does not use the notion of differentiation index, and it is applicable to any singular equation having a solution.The last three decades mathematical modeling involves problems that are modeled as systems of implicit differential equationswhereSuch systems arise in designing electric circuits (microchips for computers), mechanics (cinematics equations), nuclear energy, and others [8]. Irresolvability of the system (1) with respect to the derivativesẋ (singularity) can be a consequence of the complexity of the corresponding algebraic and transcendental transformations, as well as a consequence of the degeneration of the Jacobi matrix ∂F (ẋ, x, t)/∂ẋ on the solution {x(t)}. The statement of initial or boundary value problems for such equations should be in accordance with the system (1). The IDEs (1) are also called 'unstructured differential-algebraic equations' [1].An important particular kind of IDEs is the so-called 'structured DAEs', or simply DAEṡwhereHere initial or boundary value problems should be in accordance with the finite subsystem, and difficulties of theory and numerical solution to the initial/boundary value problems appear when the Jacobi matrix of the finite subsystem of ∂g(x, u, t)/∂u degenerates on the solution {x(t)} or {x(t), u(t)}.It is known that singular differential equations has a set of pathologies with respect to regular ODEs. So, their solutions (if they exist) can depend non-continuously on input data (initial values and parameters of equations), and the set of solutions can be infinite-dimensional (see, for example, the recent investigation [2]). Correspondingly, in the common singular case, theory and numerical methods should be created in frame of the theory of ill-posed problems [4].Often the system IDE or DAE can be transformed to the regular normal (Cauchy) by a finite number of differentiations and by algebraic transformations. The minimal number of such differentiations is called the differentiation index (DI) [8]. In applications, DAEs with indices five and higher do appear, and there are DAEs which have variable degeneracy, and they are non-reducible to the normal form.For singular differential equations with DI up to three effective specialized step-type methods exist [8]. In view of the p...