How can I make this implementation of Euler's method more efficient?

Question

I have an implementation of an Euler method for solving N-many 1st order coupled differential equations, but I feel that I did not write it as efficiently as I could, due to lack of programming experience.

Here is the implementation:

def eulerMethod(f, init0, h, om, mesh):
    """
    f          - array of ODES
    init0      - intial values     
    h          - step size
    om         - array of symbols 
    mesh       - time mesh to evolve over

    This implements Euler's method for finding numerically the
    solution of the 1st order system of N-many ODEs

    output in form of: [t:, om0:, om1:, om2:, ... ]  
    """

    numOfDE = len(f)
    t00 = mesh[0]
    soln = [[t00]]
    for i in xrange(numOfDE): # create intitial soln
        soln[0].append(init0[i])   

    subVal = {} # for diff eq substituion 
    for i in xrange(len(om)):
        subVal.update({om[i]:init0[i]})

    g = sympy.symbols('g0:%d'%len(om))
    s = sympy.symbols('s0:%d'%len(om))

    # set up dictionary for equations 
    eqDict = {g[0]:init0[0]}    
    for i in xrange(1,len(om)):
        eqDict[g[i]] = init0[i]

    for i in xrange(6): # number of steps
        for i in xrange(len(om)): # performs euler steps
            eqDict[s[i]] = eqDict[g[i]] + h*1.0*(f[i].subs(subVal))

        for i in xrange(len(om)): # set recursive values
            eqDict[g[i]] = eqDict[s[i]]

        t00 += h 
        soln.append([t00])
        for i in xrange(numOfDE): # append rest of solutions 
            soln[len(soln)-1].append(eqDict[s[i]])

        subVal = {} # update values to be subsititied        
        for i in xrange(len(om)):
            subVal.update({om[i]:eqDict[g[i]]})

    return soln       

I know my naming is sort of confusing, but I just wanted to see how solid my algorithm is. I will be using this to typically solve 1000 coupled differential equations.

Answer 1

• Naming is confusing indeed. What is sy, om, etc?

• I do not understand a need for eqDict and subVal. Using init directly is much more straightforward:

for example,

 result = [ init0[i] + h*f[i].subs(init0) for i in xrange(init0) ]
 init0 = [ x for x in result]

• Updating subVal at the end of the function is meaningless.

• It is very unclear why the Euler step is repeated 6 times for one time step.

• Is there a reason to pass mesh, if the code only uses mesh[0]?

• The description should mention that the function only addresses autonomous systems.