#  Python 3
#  File: power_demo.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

"""  
    Add up terms in a power series
  
       f(x) = sum a_n x^n

       The sum includes terms from n = 0 to n = n_max 
    Compute the coeffients a_n, put them in a numpy array,
    then add up the terms in the sum.
"""

import numpy as np    # load the nympy library, call it np

print("Demo for Differential Equations!")
print("Using a loop to add up terms in a power series")
print("Use a numpy one-index array to store the power series coefficients")

#  Compute the power series coefficients for the exponential function

n_max = 20    # take terms x^n with n up to and including n_max
x     = 4.    # the floating point value of x for the power series

a = np.zeros(n_max+1)   # coefficients for n = 0,1,...,n (n+1 in all)

#   Calculate the coefficients
n_fac = 1.   # a float, will be n! (n factorial)
for n in range(n_max+1):
    a[n] = 1./n_fac         # a_n = 1/n!
    n_fac = n_fac*(n+1)  # calculate (n+1)! from n!

#      output formatting is complicated.  The AI can do it for you
    out_string = "n = {n:3d},  a_n = {a_n:10.3e}".format(n=n, a_n = a[n])
    print(out_string)

print("The array of coefficients is " + str(a))

# compute the power series sum a_n*x*n

x_n = 1.   # will be x^n, start with x^0=1
sum = 0.   # float, will be the sum
for n in range(n_max+1):   # python range(L) goes from 0 to L-1
    sum = sum + a[n]*x_n   # add in term n in the sum
    x_n *= x            # calculate x^(n+1) 

out_string = "The sum is {sum:15.11f}, using terms up to n = {n_max:3d}"
out_string = out_string.format( sum=sum, n_max = n_max)
print(out_string)

out_string = "exp(x) is  {ex:15.11f}, and x = {x:5.2f}"
out_string = out_string.format(ex = np.exp(x), x = x)
print(out_string)