text: import numpy as np
  
  def forward_euler(f,y0,Delta_t,numsteps):
      """Perform numsteps of the forward euler method starting at y0
      of the ODE y'(t) = f(y,t)
      Args:
          f: function to integrate takes arguments y,t
          y0: initial condition
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      y = np.zeros(numsteps+1)
      t = np.arange(numsteps+1)*Delta_t
      y[0] = y0
      for n in range(1,numsteps+1):
          y[n] = y[n-1] + Delta_t * f(y[n-1], t[n-1])
      return t, y
  
  def inexact_newton(f,x0,delta = 1.0e-7, epsilon=1.0e-6, LOUD=False):
      """Find the root of the function f via Newton-Raphson method
      Args:
          f: function to find root of
          x0: initial guess
          delta: finite difference parameter
          epsilon: tolerance
          
      Returns:
          estimate of root
      """
      x = x0
      if (LOUD):
          print("x0 =",x0)
      iterations = 0
      while (np.fabs(f(x)) > epsilon):
          fx = f(x)
          fxdelta = f(x+delta)
          slope = (fxdelta - fx)/delta
          if (LOUD):
              print("x_",iterations+1,"=",x,"-",fx,"/",slope,"=",x - fx/slope)
          x = x - fx/slope
          iterations += 1
      #print("It took",iterations,"iterations")
      return x #return estimate of root
  
  def backward_euler(f,y0,Delta_t,numsteps):
      """Perform numsteps of the backward euler method starting at y0
      of the ODE y'(t) = f(y,t)
      Args:
          f: function to integrate takes arguments y,t
          y0: initial condition
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      y = np.zeros(numsteps+1)
      t = np.arange(numsteps+1)*Delta_t
      y[0] = y0
      for n in range(1,numsteps+1):
          solve_func = lambda u: u-y[n-1] - Delta_t*f(u,t[n])
          y[n] = inexact_newton(solve_func,y[n-1])
      return t, y
  
  def inexact_newton(f,x0,delta = 1.0e-7, epsilon=1.0e-6, LOUD=False):
      """Find the root of the function f via Newton-Raphson method
      Args:
          f: function to find root of
          x0: initial guess
          delta: finite difference parameter
          epsilon: tolerance
          
      Returns:
          estimate of root
      """
      x = x0
      if (LOUD):
          print("x0 =",x0)
      iterations = 0
      while (np.fabs(f(x)) > epsilon):
          fx = f(x)
          fxdelta = f(x+delta)
          slope = (fxdelta - fx)/delta
          if (LOUD):
              print("x_",iterations+1,"=",x,"-",fx,"/",slope,"=",x - fx/slope)
          x = x - fx/slope
          iterations += 1
      #print("It took",iterations,"iterations")
      return x #return estimate of root
  
  def backward_euler(f,y0,Delta_t,numsteps):
      """Perform numsteps of the backward euler method starting at y0
      of the ODE y'(t) = f(y,t)
      Args:
          f: function to integrate takes arguments y,t
          y0: initial condition
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      y = np.zeros(numsteps+1)
      t = np.arange(numsteps+1)*Delta_t
      y[0] = y0
      for n in range(1,numsteps+1):
          solve_func = lambda u: u-y[n-1] - Delta_t*f(u,t[n])
          y[n] = inexact_newton(solve_func,y[n-1])
      return t, y
  
  def RK4(f,y0,Delta_t,numsteps):
      """Perform numsteps of the 4th order Runge-Kutta method starting at y0
      of the ODE y'(t) = f(y,t)
      Args:
          f: function to integrate takes arguments y,t
          y0: initial condition
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      y = np.zeros(numsteps+1)
      t = np.arange(numsteps+1)*Delta_t
      y[0] = y0
      for n in range(1,numsteps+1):
          dy1 = Delta_t * f(y[n-1], t[n-1])
          dy2 = Delta_t * f(y[n-1] + 0.5*dy1, t[n-1] + 0.5*Delta_t)
          dy3 = Delta_t * f(y[n-1] + 0.5*dy2, t[n-1] + 0.5*Delta_t)
          dy4 = Delta_t * f(y[n-1] + dy3, t[n-1] + Delta_t)
          y[n] = y[n-1] + 1.0/6.0*(dy1 + 2.0*dy2 + 2.0*dy3 + dy4)
      return t, y
  
  def forward_euler_system(Afunc,c,y0,Delta_t,numsteps):
      """Perform numsteps of the forward euler method starting at y0
      of the ODE y'(t) = A(t) y(t) + c(t)
      Args:
          Afunc: function to compute A matrix
          c: nonlinear function of time
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      unknowns = y0.size
      y = np.zeros((unknowns,numsteps+1))
      t = np.arange(numsteps+1)*Delta_t
      y[0:unknowns,0] = y0
      for n in range(1,numsteps+1):
          yold = y[0:unknowns,n-1]
          A = Afunc(t[n-1])
          y[0:unknowns,n] = yold + Delta_t * (np.dot(A,yold) + c(t[n-1]))
      return t, y
  
  def swap_rows(A, a, b):
      """Rows two rows in a matrix, switch row a with row b
      
      args:
          A: matrix to perform row swaps on
          a: row index of matrix
          b: row index of matrix
          
      returns: nothing
      
      side effects:
      changes A to rows a and b swapped
      """
      assert (a>=0) and (b>=0)
      N = A.shape[0] #number of rows
      assert (a<N) and (b<N) #less than because 0-based indexing
      temp = A[a,:].copy()
      A[a,:] = A[b,:].copy()
      A[b,:] = temp.copy()
  def BackSub(aug_matrix,x):
      """back substitute a N by N system after Gaussian elimination
      
      Args:
          aug_matrix: augmented matrix with zeros below the diagonal
          x: length N vector to hold solution
      Returns:
          nothing
      Side Effect:
      x now contains solution
      """
      N = x.size
      for row in range(N-1,-1,-1):
          RHS = aug_matrix[row,N]
          for column in range(row+1,N):
              RHS -= x[column]*aug_matrix[row,column]
          x[row] = RHS/aug_matrix[row,row]
      return
  def GaussElimPivotSolve(A,b,LOUD=0):
      """create a Gaussian elimination with pivoting matrix for a system
      
      Args:
          A: N by N array
          b: array of length N
      Returns:
          solution vector in the original order
      """
      [Nrow, Ncol] = A.shape
      assert Nrow == Ncol
      N = Nrow
      #create augmented matrix
      aug_matrix = np.zeros((N,N+1))
      aug_matrix[0:N,0:N] = A
      aug_matrix[:,N] = b
      #augmented matrix is created
      
      #create scale factors
      s = np.zeros(N)
      count = 0
      for row in aug_matrix[:,0:N]: #don't include b
          s[count] = np.max(np.fabs(row))
          count += 1
      if LOUD:
          print("s =",s)
      if LOUD:
          print("Original Augmented Matrix is\n",aug_matrix)
      #perform elimination
      for column in range(0,N):
          
          #swap rows if needed
          largest_pos = np.argmax(np.fabs(aug_matrix[column:N,column]/s[column])) + column
          if (largest_pos != column):
              if (LOUD):
                  print("Swapping row",column,"with row",largest_pos)
                  print("Pre swap\n",aug_matrix)
              swap_rows(aug_matrix,column,largest_pos)
              #re-order s
              tmp = s[column]
              s[column] = s[largest_pos]
              s[largest_pos] = tmp
              if (LOUD):
                  print("A =\n",aug_matrix)
          #finish off the row
          for row in range(column+1,N):
              mod_row = aug_matrix[row,:]
              mod_row = mod_row - mod_row[column]/aug_matrix[column,column]*aug_matrix[column,:]
              aug_matrix[row] = mod_row
      #now back solve
      x = b.copy()
      if LOUD:
          print("Final aug_matrix is\n",aug_matrix)
      BackSub(aug_matrix,x)
      return x
  
  def backward_euler_system(Afunc,c,y0,Delta_t,numsteps):
      """Perform numsteps of the forward euler method starting at y0
      of the ODE y'(t) = A(t) y(t) + c(t)
      Args:
          Afunc: function to compute A matrix
          c: nonlinear function of time
          y0: initial condition
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      unknowns = y0.size
      y = np.zeros((unknowns,numsteps+1))
      t = np.arange(numsteps+1)*Delta_t
      y[0:unknowns,0] = y0
      for n in range(1,numsteps+1):
          yold = y[0:unknowns,n-1]
          A = Afunc(t[n])
          LHS = np.identity(unknowns) - Delta_t * A
          RHS = yold + c(t[n])*Delta_t
          y[0:unknowns,n] = GaussElimPivotSolve(LHS,RHS)
      return t, y
  
  def cn_system(Afunc,c,y0,Delta_t,numsteps):
      """Perform numsteps of the forward euler method starting at y0
      of the ODE y'(t) = A(t) y(t) + c(t)
      Args:
          Afunc: function to compute A matrix
          c: nonlinear function of time
          y0: initial condition
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      unknowns = y0.size
      y = np.zeros((unknowns,numsteps+1))
      t = np.arange(numsteps+1)*Delta_t
      y[0:unknowns,0] = y0
      for n in range(1,numsteps+1):
          yold = y[0:unknowns,n-1]
          A = Afunc(t[n])
          LHS = np.identity(unknowns) - 0.5*Delta_t * A
          A = Afunc(t[n-1])
          RHS = yold + 0.5*Delta_t * np.dot(A,yold) + 0.5*(c(t[n-1]) + c(t[n]))*Delta_t
          y[0:unknowns,n] = GaussElimPivotSolve(LHS,RHS)
      return t, y
  
  def RK4_system(Afunc,c,y0,Delta_t,numsteps):
      """Perform numsteps of the forward euler method starting at y0
      of the ODE y'(t) = f(y,t)
      Args:
          f: function to integrate takes arguments y,t
          y0: initial condition
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      unknowns = y0.size
      y = np.zeros((unknowns,numsteps+1))
      t = np.arange(numsteps+1)*Delta_t
      y[0:unknowns,0] = y0
      for n in range(1,numsteps+1):
          yold = y[0:unknowns,n-1]
          A = Afunc(t[n-1])
          dy1 = Delta_t * (np.dot(A,yold) + c(t[n-1])) 
          A = Afunc(t[n-1] + 0.5*Delta_t)
          dy2 = Delta_t * (np.dot(A,y[0:unknowns,n-1] + 0.5*dy1) 
                           + c(t[n-1] + 0.5*Delta_t))
          dy3 = Delta_t * (np.dot(A,y[0:unknowns,n-1] + 0.5*dy2) 
                           + c(t[n-1] + 0.5*Delta_t))
          A = Afunc(t[n] + Delta_t)
          dy4 = Delta_t * (np.dot(A,y[0:unknowns,n-1] + dy3) + c(t[n]))
          y[0:unknowns,n] = y[0:unknowns,n-1] + 1.0/6.0*(dy1 + 2.0*dy2 + 2.0*dy3 + dy4)
      return t, y
  
  def new2(f,y0,Delta_t,numsteps):
      """Perform numsteps of the backward euler method starting at y0
      of the ODE y'(t) = f(y,t)
      Args:
          f: function to integrate takes arguments y,t
          y0: initial condition
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      y = np.zeros(numsteps+1)
      t = np.arange(numsteps+1)*Delta_t
      y[0] = y0
      for n in range(1,numsteps+1):
          tcur = t[n-1]
          fshift = lambda u,t: f(u,t+tcur)
          tmp,y1 = backward_euler(fshift,y[n-1],Delta_t*0.5,1)
          tmp,y2 = backward_euler(fshift,y[n-1],Delta_t,1)
          y[n] = y[n-1] + Delta_t*(1.0/(1.0+Delta_t) * f(y1[1],t[n-1]+0.5*Delta_t) + 
                                   Delta_t/(1.0+Delta_t) * f(y2[1],t[n]))
      return t, y
  
  def new2_system(Afunc,c,y0,Delta_t,numsteps):
      """Perform numsteps of the forward euler method starting at y0
      of the ODE y'(t) = A(t) y(t) + c(t)
      Args:
          Afunc: function to compute A matrix
          c: nonlinear function of time
          y0: initial condition
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      unknowns = y0.size
      y = np.zeros((unknowns,numsteps+1))
      t = np.arange(numsteps+1)*Delta_t
      y[0:unknowns,0] = y0
      for n in range(1,numsteps+1):
          tcur = t[n-1]
          Af = lambda t: Afunc(t+tcur)
          cf = lambda t: c(t+tcur)
          tmp,y1 = backward_euler_system(Af,cf,y[0:unknowns,n-1],Delta_t*0.5,1)
          tmp,y2 = backward_euler_system(Af,cf,y[0:unknowns,n-1],Delta_t,1)
          y1 = y1[0:unknowns,1]
          y2 = y2[0:unknowns,1]
          th = t[n-1] + 0.5*Delta_t
          Ah = Afunc(th)
          A = Afunc(t[n])
          y[0:unknowns,n] = y[0:unknowns,n-1] + Delta_t*(1.0/(1.0+Delta_t) * (np.dot(Ah,y1) + c(th)) + 
                                   Delta_t/(1.0+Delta_t) * (np.dot(A,y2) + c(t[n])) )
      return t, y
  
  def crank_nicolson(f,y0,Delta_t,numsteps):
      """Perform numsteps of the backward euler method starting at y0
      of the ODE y'(t) = f(y,t)
      Args:
          f: function to integrate takes arguments y,t
          y0: initial condition
          Delta_t: time step size
          numsteps: number of time steps
          
      Returns:
          a numpy array of the times and a numpy
          array of the solution at those times
      """
      numsteps = int(numsteps)
      y = np.zeros(numsteps+1)
      t = np.arange(numsteps+1)*Delta_t
      y[0] = y0
      for n in range(1,numsteps+1):
          solve_func = lambda u: u-y[n-1] - 0.5*Delta_t*(f(u,t[n])
                                                         + f(y[n-1],t[n-1]))
          y[n] = inexact_newton(solve_func,y[n-1])
      return t, y