solving ordinary differential equations if odespy_methods is not listbut just the name of single odespy solverwe wrap that name in list so we always have odespy_methods as list if type(odespy_methods!type([])odespy_methods [odespy_methodsmake list of solver objects solvers [method(ff_args=[omega]for method in odespy_methodsfor solver in solverssolver set_initial_condition([ x_ ]compute the time points where we want the solution dt float(dtavoid integer division n_t int(round( /dt)time_points linspace( n_t*dtn_t+ legends [for solver in solverssolt solver solve(time_pointsv sol[:, sol[:, plot only the last periods *time_intervals_per_period no time steps to plot plot( [- :] [- :]hold('on'legends append(solver name()xlabel(' 'plot exact solution too plot( [- :]x_ *cos(omega* )[- :]' --'legends append('exact'legend(legendsloc='lower left'axis([ [- ] [- ]- *x_ *x_ ]title('simulation of % periods with % intervals per period(number_of_periodstime_intervals_per_period)savefig('tmp pdf')savefig('tmp png'show( new feature in this code is the ability to plot only the last periodswhich allows us to perform long time simulations and watch the end results without cluttered plot with too many periods the syntax [- :plots the last elements in ( negative index in python arrays/lists counts from the endwe may compare heun' method (or equivalently the rk methodwith the euler-cromer schemecompare(odespy_methods=[odespy heunodespy eulercromer]omega= x_ = number_of_periods= time_intervals_per_period= figure shows how heun' method (the blue line with small diskshas considerable error in both amplitude and phase already after - periods (upper left)but using three times as many time steps makes the curves almost equal (upper