BeginPackage["MidPoint`"] MidPoint::usage="MidPoint[{e1,e2,...},{y1,y2,...},{a1,a2,...},{t1,dt}] numerically integrates the ei as functions of the yi with initial values ai. The integration proceeds in steps of dt from 0 to t1 using the two-step MidPoint method." Begin["`Private`"] MidPoint[f_List, y_List, y0_List, {t1_, dt_}] := Module[{}, midptmeth = Append[{y0}, y0 + dt N[ f/. Thread[ y-> y0] ]]; Do[midptmeth = Append[midptmeth, midptmeth[[i-2]] + 2*dt*N[ f /. Thread[y->midptmeth[[i-1]] ]]], {i,3,Round[t1/dt]+1}]; Return[midptmeth] ] End[] EndPackage[]