Physiology models written in a description language such as CellML are becoming a popular method to handle complex cellular physiological models in biological function simulations. However, in order to fully simulate a model, boundary conditions and ordinary differential equation (ODE) solving schemes have to be combined with it. Though the former can be described in CellML, it is difficult to explicitly specify ODE solving schemes using existing tools. In this study, we defined an ODE solving scheme description language based on XML and proposed a code generation system for biological function simulations. By using the proposed system, biological simulation programs using various ODE solving schemes can be easily generated. We designed a two-stage approach where the system generates a set of equation associating with the physiological model variable values at a certain time t with values at t plus delta t in the first stage and generates the programs calculating the time evolution of the model in the second stage. This approach enables the flexible construction of code generation modules that can support complex sets of formulas. We evaluated the relationship between models and their calculation accuracies by simulating complex biological models using various ODE solving schemes. Using the FHN model simulation, results showed good qualitative and quantitative correspondence with the theoretical predictions. Results for the LuoRudy1991 model showed that only first order precision was achieved. In addition, running the generated code in parallel on GPU made it possible to speed up the calculations by a factor of 50.