# -*- coding: utf-8 -*-
"""
Created on Mon Mar  8 12:12:15 2021

@author: turinici
"""

import numpy as np
import matplotlib.pyplot as plt
#%matplotlib inline
%matplotlib auto

# np.sqrt(2.)*np.sqrt(2.) == 2.
# False !!!!

#parameters
T=10.0 # final time
N=100 # number of time steps
h = T/N # 
U0 = 0.0

def ftxt(t,x):
# define the function entering the ODE x'(t) = f(t,x(t))
    return(2.0*np.sqrt(np.abs(x)))

def Explicit_Euler(h,N,ode_function,initial_value):
    """
    Parameters
    ----------
    h : time step
    N : number of time steps (integer)
    input_function : function to integrate
    initial_value : initial value of the ODE

    Returns
    -------
    List of approximate solution obtained by Explicit Euler scheme of step h.

    """

    U= [0.0]*(N+1) # creates a list on N+1 elements, filled with value 0
    U[0]=initial_value

    for ii in range(N):
        U[ii+1]=U[ii]+h*ode_function(ii*h,U[ii])
    return(U)
    
# define the time grid    
trange = np.linspace(0,T,num=N+1,endpoint=True)   
# solve by EE the ODE : use EE at any time step and put in a numpy array : solution
solution = Explicit_Euler(h,N,ftxt,U0)
solution2 = Explicit_Euler(h,N,ftxt,np.sqrt(2.)*np.sqrt(2.)-2.)

# plot the results
plt.figure(1)
#plt.plot(solution)
plt.plot(trange,solution,'*r',trange,solution2,'ob')
plt.xlabel('time (t)')
plt.ylabel('x(t)')
plt.legend(['infinite precision solution x(t)','finite precision solution x(t)'])