# Author: Alex Gezerlis

# Numerical Methods in Physics with Python (CUP, 2020)

from gauelim_pivot import gauelim_pivot

import numpy as np

def matsetup(a,b,n):

h = (b-a)/(n-1)

xs = a + np.arange(n)*h

A = np.zeros((n,n))

np.fill_diagonal(A, -2 + 30*h**2/(1-xs**2))

A[0,0] = 1; A[-1,-1] = 1

np.fill_diagonal(A[1:,:], 1 + h*xs[1:]/(1-xs[1:]**2))

A[-1,-2] = 0

np.fill_diagonal(A[:,1:], 1 - h*xs/(1-xs**2))

A[0,1] = 0

bs = np.zeros(n)

bs[0] = 0.0926587109375

bs[-1] = 0.11177050858750004

return A, bs

def riccati(a,b,n):

A, bs = matsetup(a, b, n)

ws = gauelim_pivot(A, bs)

return ws

if __name__ == '__main__':

a, b, n = 0.05, 0.49, 400

ws = riccati(a, b, n); print(ws)