title: Solving ODE in complex domain with Python (or Matlab)
body: <p>As a test for a more complicated system, I want to solve a differential equation dw/dz = w where the function w = w(z) is complex valued and z = x+iy as usual. The boundary conditions are w = i when z = i. The solution is of course complex and defined on the argand plane. I was hoping to solve this with some standard ODE solvers in python. My method is to first define a grid in the argand plane (lines of constant x and y) and then loop through each grid line and call an ODE solver on each iteration. In the below code I am attempting to integrate my differential equation between 1j and 2j, but the resulting vector of w is just 1j! Can anyone advise me what to do? Thanks</p>
  
  <pre><code>    from scipy.integrate import ode
      import numpy as np
      from matplotlib.pylab import *
  
  
      def myodeint(func, w0, z):
          w0 = np.array(w0, complex)
          func2 = lambda z, w: func(w, z)   # odeint has these the other way :/
          z0 = z[0]
          solver = ode(func2).set_integrator('zvode').set_initial_value(w0, z0)
          w = [solver.integrate(zp) for zp in z[1:]]
          w.insert(0, w0)
          return np.array(w)
  
      def func2(w, z, alpha):
          return alpha*w
  
      if __name__ == '__main__':
          # Set grid size in z plane
          x_max = 3
          x_min = 0
          y_max = 3
          y_min = 0
  
          # Set grid resolution
          dx = 0.1
          dy = 0.1
  
          # Number of nodes
          x_nodes = int(np.floor((x_max-x_min)/dx)+1)
          y_nodes = int(np.floor((y_max-y_min)/dy)+1)
  
          # Create array to store value of w(z) at each node
  
          ww = np.zeros((y_nodes,x_nodes), complex)
  
          # Set boundary condition: w = w0 at x = x0, y = y0
          x0 = 0
          y0 = 1
          i0 = (x0-x_min)/dx
          j0 = (y_max-y0)/dy
          w0 = 1j
          ww[j0,i0] = w0
  
  
          z0 = 1j
  
          alpha = 1
          z = np.linspace(z0, z0+1j, 200)
          w = myodeint(lambda w, z: func2(w, z, alpha), [w0, 0, 0], z)
  </code></pre>