#  Python 3
#  File: Euler_saddle.py
#  Demo Python code for Differential Equations class, Spring 2026
#  Authors: Jonathan Goodman, goodman@cims.nyu.edu
#  class web site: https://math.nyu.edu/~goodman/teaching/DifferentialEquations2026/DifferentialEquations.html

"""  
    Solve two component ODE systems and make trajectory plots in the
    phase plane.  Use the forward Euler solver
    This example, a linear ODE system that makes trajectories for a saddle point
 """

import numpy             as np    # load the nympy library, call it np
import matplotlib.pyplot as plt   # the plotting package

print("Euler method for trajectories in a phase plane.")

#  Set parameters

x_0 = np.array([1.,0.])    # initial condition
T   = 4.    # integrate to this time
A   = np.array([[-.7 , 2.],
                [ .4, 1.6]])

n = 10000  # Number of time steps for each trajectory

x_0_list = [ np.array([ 1.  , 1.]) ,      # a Python list of initial values
             np.array([-1.  , 1.]) ,
             np.array([ 1.  , -.2]),
             np.array([-1.  , .1]) ,
             np.array([-.7  , 0.6]) ,
             np.array([-.988, 0.152]) ,
          ]
d = len(x_0_list[0])         # d = problem dimension, found from any x_0

n_runs    = len(x_0_list)    # the number trajectories computed

#  Use the integrated eigenvalue solver to find the eigenvalues
#  and eigenvectors of A

eigenvalues, eigenvectors = np.linalg.eig(A)
for k in range(d):
   lam = eigenvalues[k]          # not lambda because it's a keyword in Python
   v      = eigenvectors[:,k]    # column k of the eivenvector matrix
   eig_info = "lambda = {lam:11.4e} has eigenvector ({v_1:11.4e},{v_2:11.4e})"
   eig_info = eig_info.format( lam = lam, v_1 = v[0], v_2=v[1])
   print(eig_info)

def f(x,A):
   """
      The ODE is (dx/dt) = Ax
      Inputs:
         x (numpy array of floats): the state, d components
         A (numpy 2 index square array: the d by d matrix
      Return value
         (numpy float array): computed value of f
   """
   return A@x    #   @ is matrix/vector or matrix/matrix multiplication

fig,ax = plt.subplots()     # create the figure, empty at first, where plots will go
ax.set_title(r"Euler solver, circles, to $T=${T:4.2f}".format(T=T))

#    Compute the function curves one by one

i = 0                              # index that determines which n is being used
for x_0 in x_0_list:                  # take each n value in turn

   dt = T/n                        # n time steps of size dt go from 0 to T
   
   t_vals = np.zeros(n+1)          # hold t_0=0 up to T = t_n
   x_vals = np.zeros([2,n+1])      # hold x_n from x_0 (initial val) to x_n
                                   # a 2Xn array, one column for each x value
   
   t = 0.                          # starting time
   x = x_0                         # initial condition
   
   t_vals[0] = t                   # not necessary because it's already zero
   x_vals[:,0] = x_0
   ax.plot(x_0[0],x_0[1], "o")     # put a big dot at the beginning of the curve
   
   for k in range(n):              # the time step loop
      x = x + dt*f(x,A)            # the ODE, move from x_k to x_{k+1}
      t = t + dt                   # replace t_k with t_{k+1}
      
      t_vals[k+1]   = t            # these are the values at the next time
      x_vals[:,k+1] = x
   
   i           = i+1
      
   curve_label = "curve = {number:2d}".format(number = i)
   x_1 = x_vals[0,:]
   x_2 = x_vals[1,:]
   ax.scatter( x_1, x_2, label = curve_label, s = .1)   # line in the v_1 direction

plot_scale = 2.
plt.gca().set_aspect('equal')

#   Put the eigenvector directions into the plot

v_1 = eigenvectors[:,0]
e_s = np.array([-2*v_1[0],2*v_1[0]])
e_e = np.array([-2*v_1[1],2*v_1[1]])
ax.plot( e_s, e_e)

v_2 = eigenvectors[:,1]
e_s = np.array([-2*v_2[0],2*v_2[0]])
e_e = np.array([-2*v_2[1],2*v_2[1]])
ax.plot( e_s, e_e)

ax.set_xlim( -plot_scale, plot_scale)     # make the curves go into a square
ax.set_ylim( -plot_scale, plot_scale)
ax.grid()
ax.set_xlabel(r"$x_1$")
ax.set_ylabel(r"$x_2$")
plt.savefig("Euler_saddles.pdf")