# A script to illustrate the backward/implicit Euler method for systems or scalar ODEs.
# Key points:
# The spatial variables are y[0],...,y[d-1] for a system with d dependent variables.
# t0 is the initial time, T is the final time
# The vector "f" is an expression involving t,y[0],...,y[d-1] that defines y'=f(t,y)
# Define f as below (scalar and vector examples shown).
# Also Jf is the Jacobian matrix of f with respect to y, supplied by user.
# Solution times are returned in row vector "tv"
# Solution values are return in matrix "yv". Rows correspond to times, columns
# to variables y[0],...,y[d-1].

reset();
var('t y');

# SYSTEM EXAMPLE: Define f in y' = f(t,y). Initial condition is y(t0) = ic.
# You must also define "Jf", the Jacobian matrix for f with respect to y (used in Newton's
# method for solving the implicit equation).
def f(t,y):
    return vector([-sin(t)*y[1],-y[0]+cos(y[1])])
def Jf(t,y):
    return(matrix([[0,-sin(t)],[-1,-sin(y[1])]]))
ic = vector(SR,[1.0,2.0]);

#SCALAR EXAMPLE. Uncomment to test.
#def f(t,y):
#    return vector([-sin(y[0])])
#def Jf(t,y):
#    return(matrix([[-cos(y[0])]]))
#ic = vector(SR,[1.0]);

# Set solution interval t0 <= t <= T, and stepsize.
t0 = 0.0;
T = 1.0;
stepsize = 0.1;

#Call routine for implicit Euler's method.
load('implicit_euler_sys.sage');
[tv4,yv4] = implicit_euler_sys(f,Jf,ic,t0,T,stepsize)

#Print solution values at final time.
N = len(tv4)-1;
d = len(ic);
print('Implicit Euler Method');
print('time',tv4[N]);
for j in range(d):
    print('y[',j,'] = ',yv4[N,j]);

#Plot first component of solution versus time, as simple example.
yp = yv4.column(0);
p = line(list(zip(tv4,yp)),rgbcolor=[1,0,1],legend_label='Implicit Euler')
p.set_legend_options(loc='upper right');
p.axes_labels(['$t$','$y0(t)$'])
show(p)