An interesting property of the inverse F-transform f ^ of a continuous function f on a given interval [ a , b ] says that the integrals of f ^ and f on [ a , b ] coincide. Furthermore, the same property can be established for the restrictions of the functions to all subintervals [ a , p k ] of the fuzzy partition of [ a , b ] used to define the F-transform. Based on this fact, we propose a new method for the numerical solution of ordinary differential equations (initial-value ordinary differential equation (ODE)) obtained by approximating the derivative x · ( t ) via F-transform, then computing (an approximation of) the solution x ( t ) by exact integration. For an ODE, a global second-order approximation is obtained. A similar construction is then applied to interval-valued and (level-wise) fuzzy differential equations in the setting of generalized differentiability (gH-derivative). Properties of the new method are analyzed and a computational section illustrates the performance of the obtained procedures, in comparison with well-known efficient algorithms.