# solves mixed problem with -u''(x)=f(x), 0<x<1, u(0)=0, u'(1)=sigma
# user specified functions u(x) and grid x as input to code

import numpy as np
import matplotlib.pyplot as plt

# 
def f(x):
    return (np.pi**2)*np.sin(np.pi*x)

def U_neu(x,sigma):
    n = x.shape[0]
    h = np.diff(x)
    A = np.zeros((n,n))
    b = np.zeros((n,1))
    # loop over elements
    for i in range(n-1):
        A[i:i+2][:,i:i+2] = A[i:i+2][:,i:i+2] + np.array([[1,-1],[-1,1]])/h[i]
        # use midpoint ('barycentric') one node quadrature
        b[i:i+2] = b[i:i+2] + f(np.mean(x[i:i+2]))*h[i]/2
    # implement boundary conditions
    A[0] = np.zeros((1,n))
    A[0][0] = 1
    b[0] = 0
    b[-1] = b[-1] + sigma
    # solve system
    u = np.linalg.solve(A,b)
    return u

sigma=-np.pi
x = np.linspace(0,1,100)
u=U_neu(x,sigma)   
plt.plot(x,u)