title: Euler Method implementation in Python gives a stable result but it should be unstable
body: <p>I am trying to solve this differential equation with the Euler Method using Python3:</p>
  
  <p><a href="https://i.stack.imgur.com/QzkhH.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/QzkhH.png" alt="enter image description here"></a></p>
  
  <p>According to Wolfram Alpha, that's the plot of the correct equation. </p>
  
  <p><a href="https://i.stack.imgur.com/DsbU9.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/DsbU9.png" alt="enter image description here"></a></p>
  
  <p>Again, according to Wolfram Alpha, in this case, the classic Euler Method should <strong>not</strong> be stable, as you can see by the end of the interval:</p>
  
  <p><a href="https://i.stack.imgur.com/ELvBB.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/ELvBB.png" alt="enter image description here"></a></p>
  
  <p>However, on my implementation, Euler method provides a stable result, which <strong>is strange</strong>. I wonder that my implementation is wrong for some reason. Nonetheless, I can't find the error.</p>
  
  <p>I generated some points and a plot comparing my approximation and the analytic output of the function. In blue, the analytic result as control group. In red, the output of my implementation:</p>
  
  <p><a href="https://i.stack.imgur.com/IZPVT.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/IZPVT.png" alt="enter image description here"></a></p>
  
  <p>That's my code:</p>
  
  <pre><code>import math
  import numpy as np
  from matplotlib import pyplot as plt
  import pylab
  
  def f(x):
  
      return (math.e)**(-10*x)
  
  def euler(x):
  
      y_init = 1
      x_init = 0
  
      old_dy_dx = -10*y_init
  
      old_y = y_init 
  
      new_y = None
  
      new_dy_dx = None
  
      delta_x = 0.001
  
      limite = 0
  
      while x&gt;limite:
  
          #for i in range(1,6):
  
          new_y = delta_x*old_dy_dx + old_y
          #print ("new_y", new_y)
  
          new_dy_dx = -10*new_y
          #print ("new dy_dx", new_dy_dx)
  
          old_y = new_y
          #print ("old_y", old_y)
  
          old_dy_dx = new_dy_dx
          #print ("old delta y_delta x", old_dy_dx)
          #print ("iterada",i)
  
          limite = limite +delta_x
  
      return new_y
  
  t = np.linspace(-1,5, 80)
  
  lista_outputs = []
  
  for i in t:
      lista_outputs.append(euler(i))
      print (i)
  
  # red dashes, blue squares and green triangles
  plt.plot(t, f(t), 'b-', label='Output resultado analítico')
  plt.plot(t , lista_outputs, 'ro', label="Output resultado numérico")
  plt.title('Comparação Euler/Analítico - tolerância: 0.3')
  pylab.legend(loc='upper left')
  plt.show()
  </code></pre>
  
  <p>Thanks for the help.</p>
  
  <p>============================================================</p>
  
  <p><strong>UPDATE</strong></p>
  
  <p>With @SourabhBhat help, I was able to see that my implementation was, actually, right. It was, indeed, generating an instability. Besides increasing the step size, I needed to do some zoom in to see it happening.</p>
  
  <p>The picture bellow speaks for itself (step size of 0.22):</p>
  
  <p><a href="https://i.stack.imgur.com/na1FS.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/na1FS.png" alt="enter image description here"></a></p>