function [t,y] = bvpsolve(u,v,w,a,b,alpha,beta,m)

h = (b-a)/(m+1);
t = linspace(a,b,m+2);

% Set up subdiagonal, diagonal, and superdiagonal of mxm system matrix
subdiag(1:m-1) = -1 - feval(w,t(3:m+1))*h/2;
d(1:m) = 2 + h^2*feval(v,t(2:m+1));
superdiag(1:m-1) = -1 + feval(w,t(2:m))*h/2;

% Set up right-hand side
f(1) = -h^2*feval(u,t(2)) - alpha*(-1 - feval(w,t(2))*h/2);
f(2:m-1) = -h^2*feval(u,t(3:m));
f(m) = -h^2*feval(u,t(m+1)) - beta*(-1 + feval(w,t(m+1))*h/2);

% Solve tridiagonal system
y = zeros(m+2,1);
y(1) = alpha;
y(m+2) = beta;

y(2:m+1) = tridiag(subdiag,d,superdiag,f);

