%tpbvp.m -- Solve boundary value problem
% y''=py'+qy+r
% y(a)=g1, y(b)=g2
clear all
N=20
a=0
b=1
g1=0
g2=log(2)
p=@(x) -2*x./(1+x.^2)
q=@(x) 1+0*x
r=@(x) 2./(1+x.^2)-log(1+x.^2)
h=(b-a)/N
xn=[a:h:b]'
Aa=-(1+h/2*p(xn(3:N)))
Ab=2+h^2*q(xn(2:N))
Ac=(-1+h/2*p(xn(2:N-1)))
f=-h^2*r(xn(2:N))
f(1)=f(1)+(1+h/2*p(xn(2)))*g1
f(N-1)=f(N-1)+(1-h/2*p(xn(N)))*g2
ym=tridiag(Aa,Ab,Ac,f)
yn=[g1; ym; g2]
exact=@(x) log(1+x.^2)
yexact=exact(xn)
norm(yn-yexact)


