function [U,err]=bvpfd(J,f,uexact)
% BVPFD  Function which solves a certain type of linear ODE boundary value 
% problem by finite differences on the interval 0 <= x <= 1:
%   u'' + sin(x) u' - 5 u = f(x),  u(0)=1,  u'(1)=0.
% Format:
%   [U,err]=bvpfd(J,f,uexact)
% where
%        U = approximate solution (J+1 length row vector)
%      err = maximum error (if computable; otherwise NaN)
%        J = number of subintervals
%        f = user-supplied function handle (or inline) for right side
%  [uexact = OPTIONAL;  user-supplied fcn handle (or inline) for exact soln]    
% Example I; to solve problem 3 on A#9:
%   >> f=inline('x.^2-2','x')
%   >> U=bvpfd(20,f);
% Example II; to solve a problem with manufactured soln, to see convergence:
%   >> f=inline('-(pi^2+5)*cos(pi*x)-pi*sin(x).*sin(pi*x)','x')
%   >> uexact=inline('cos(pi*x)','x')
%   >> J=[10 20 40 80 160 250 500 1000 2000 5000];
%   >> for k=1:length(J), J(k), [U,err(k)]=bvpfd(J(k),f,uexact); err(k), end
%   >> loglog(1./J,err,'o'), xlabel('\Delta x'), ylabel('error')
% ELB 4/20/05

dx=1/J;   x=0:dx:1;

% build A by diagonal first, then superdiagonal, then subdiagonal
A=sparse(1:J,1:J,repmat(-2/(dx^2)-5,1,J),J,J);
for k=1:J-1,  A(k,k+1)=1/(dx^2)+sin(x(k+1))/(2*dx);  end
for k=2:J-1,  A(k,k-1)=1/(dx^2)-sin(x(k+1))/(2*dx);  end
A(J,J-1)=2/(dx^2);
% figure, spy(A)  % its tridiagonal

% build right hand side
b=feval(f,x(2:J+1))';  
b(1)=b(1)-1/(dx^2)+sin(x(2))/(2*dx);

% solve and display
U=A\b;
U=[1 U'];  % append left bdry cond
plot(x,U), xlabel x, ylabel u
title('approx soln to BVP  u'''' + (sin x) u'' - 5 u = f(x)')

if nargin==3
   Uex=feval(uexact,x);  err=max(abs(U-Uex));
elseif nargin == 2, err=NaN;
else, error('number of arguments should be 2 or 3'), end