  function Nov18(n)
% Illustrates the function TwoPtBVP by apply it to the problem
%         du/dx^2 + u = 2x*sin(x) - 2*cos(x),   u(0) = 0, u'(pi) = -pi
close all
a = 0;  ua = 0;
b = pi; ub = -pi;
u = TwoPtBVPMod(a,ua,b,ub,n,@p,@q,@r);
x = linspace(a,b,n)';
uExact = x.*sin(x);
err = max(abs(u-uExact));
plot(x,u,'r',x,uExact,'k')
title(sprintf('h = %10.2e    err = %10.2e',(b-a)/(n-1),err),'fontsize',15)
axis([0 3.5 0 2])
shg

  function u = TwoPtBVPMod(a,ua,b,ub,n,p,q,r)
% Assume a<b and that n > 2. Let h = (b-a)/(n-1) and define x(i) = a+(i-1)h
% for i=1:n.
% u is a column n-vector with the property that u(i) approximates the
% solution to the following BVP at x(i):
%
%   -d/dx(p(x)du/dx) + q(x)u(x) = r(x)     u'(a) = ua, u'(b) = ub;
%

% Assume 2<=i<=n-1 and discretize at x(i):
%
%    u(i+1)-u(i)            u(i)-u(i-1)
%   ------------ p(i)   -  ----------- p(i-1)  + q(i)u(i)  = r(i)
%         h                       h
% - ----------------------------------------
%                       h
% where 
%                p(i) = p(x(i)+h/2) i=1:n-1
%                q(i) = q(x(i)),    i=2:n
%                r(i) = r(x(i)),    i=2:n
% 
%  Rearranging we get the following for i=2:n-1...
%
%       -u(i-1)p(i-1) + u(i)*(p(i-1)+p(i)+h^2*q(i)) - u(i+1)p(i)  = h^2r(i)
%
%  We need to evaluate p, q, and r...

x = linspace(a,b,n);
h = (b-a)/(n-1);

pVec = p(x(1:n-1)+h/2);
qVec = q(x);
rVec = r(x);

% Noting that u(1) = ua and u(n) = ub we assemble these equations in to 
% a system Tu(2:n) = rhs where

T = zeros(n,n);
rhs = zeros(n,1);
% Left endpoint condition: (u(2)-u(1))/h = ua
T(1,1)=-1; T(1,2)=1; rhs(1) = h*ua;
for i=2:n-1 
    % The ith equation from above is
    %   -u(i-1)p(i-1) + u(i)*(p(i-1)+p(i)+h^2*q(i)) - u(i+1)p(i)  = h^2r(i)
    T(i,i) = pVec(i-1)+pVec(i)+h^2*qVec(i);
    T(i,i-1) = -pVec(i-1);
    T(i,i+1) = -pVec(i);
    rhs(i) = h^2*rVec(i);
end
% Right endpoint condition: (u(n)-u(n-1))/h = ub
T(n,n) = 1; T(n,n-1) = -1; rhs(n) = h*ub;
u = T\rhs;

function y = p(x)
y = ones(length(x),1);

function y = q(x)
y = ones(length(x),1);


function y = r(x)
y = 2*x.*sin(x) - 2*cos(x);


        
      
        



