#!/usr/bin/env python
# coding: utf-8

# In[9]:


#solve BVP example using finite differences
K = lambda x: 5 + x
Kp = lambda x: 1

#Ky'' + K'y' = 0
#y'' = -K'y'/K

a = 0
b = 10
ya = 30
yb = 20


# In[10]:


N = 5
h = (b - a) / N #grid spacing

from numpy import zeros

A = zeros((N+1, N+1))
c = zeros(N+1)


A[0, 0] = 1
c[0] = ya
A[N, N] = 1
c[N] = yb

for i in range(1, N):
    x = a + i*h
    Kx = K(x)
    Kpx = Kp(x)
    A[i, i-1] = (Kx / (h**2)) - (Kpx / (2*h))
    A[i, i] = - 2*Kx / (h**2)
    A[i, i+1] = (Kx / (h**2)) + (Kpx / (2*h))


# In[3]:


A


# In[11]:


from numpy.linalg import solve
y = solve(A, c)


# In[12]:


y


# In[13]:


from matplotlib.pyplot import plot, figure, grid, title, xlabel, ylabel, legend, show
from numpy import linspace
x = linspace(a, b, N+1)
plot(x, y, 's')


# In[14]:


x, y


# In[17]:


#solve using the numerical solver solve_bvp
from scipy.integrate import solve_bvp
from numpy import array

#new variables for first-order system: z[0] = y, z[1]= y'
f = lambda x, z: [z[1], -(Kp(x)/K(x))*z[1]] #differential equation(s)
bc = lambda za, zb: [za[0] - ya, zb[0] - yb] #boundary conditions
x = [a, b] #domain
s = (yb-ya)/(b-a) #average slope in y between a and b
z = [[ya, yb], [s, s]] #known or guessed values at a and b 
                       #for z[0]=y and z[1]=y'
sol = solve_bvp(f, bc, x, z)


# In[18]:


sol


# In[20]:


plot(sol.x, sol.y[0, ], 's')


# In[21]:


x_plot = linspace(a, b, 100) #denser grid of points
y_plot = sol.sol(x_plot)[0]
plot(x_plot, y_plot)


# In[ ]: