function [ERR, YLEFT] = bvp3diff(N)
%BVP3DIFF  Solve boundary value problem 3 in section 2.4
% Try "[err yl]=bvp3diff(5);".  
% ERR=maximum of the error by comparing to the exact solution.
% YLEFT=approximation of y(-1).
% (Ed Bueler; 2/4/02)

h = 1/N;  x=-1+h:h:1-h; % note we want x_0=-1 and x_1=-1-h
A=zeros(2*N-1,2*N-1);  b=zeros(1,2*N-1);

A(1,1:3)=[-4-9*h*x(1), 2+10*h*x(1), 2-h*x(1)];
b(1)=11*h*h*x(1);
for j=2:2*N-2
   A(j,j-1:j+1)=[1-.5*h*x(j), -2, 1+.5*h*x(j)];
   b(j)=h*h*x(j);
end;
A(2*N-1,2*N-2:2*N-1)=[1-.5*h*x(2*N-1), -2];
b(2*N-1)=h*x(2*N-1)*(h-.5)-1;

figure(2);spy(A);
title('Note: not quite tridiagonal!');figure(1);

y=(A\b')'; % computed at j=1,...,2*N-1 only

yexact=x+exp(.5)*sqrt(pi/2)*(erf(1/sqrt(2))-erf(x/sqrt(2)));
ERR=max(abs(y-yexact));

XP=-1:h:1;
YP(2:2*N)=y; YP(2*N+1)=1;
YP(1)=(18*YP(2)-9*YP(3)+2*YP(4))/11; % use (2.36) to get y(x_0)
YLEFT=YP(1); % return y-value of free end
figure(1); plot(XP,YP,':'); hold;

XEP=-1:.01:1;
YEXACTP=XEP+exp(.5)*sqrt(pi/2)*(erf(1/sqrt(2))-erf(XEP/sqrt(2)));
plot(XEP,YEXACTP); 
title('Approximate (dashed) and exact (solid) solution of BVP 3.');
hold off;