% finite difference approximation of the solution to the BVP
% u''(x) + (1/x)u'(x) + u(x) = x^2,  a < x < b,
%                       u(a) = 0,
%                       u(b) = 0.

a = 10; b = 20;    % interval [a,b]
N = 100;            
h = (b-a)/(N+1);

D1 = -diag(ones(N,1)) + diag(ones(N-1,1),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% The above is the derivative matrix (times h) 
%%% we derived in class. Here's a better one:
%D1 = .5*(diag(ones(N-1,1),1) - diag(ones(N-1,1),-1) );
%%% This comes from a symmetric difference quotient of the form
%%%  u'(x_j) (approx)= [u(x_{j+1}) - u(x_{j-1})]/(2h),
%%% which has error of order h^2 instead of the order h we get with
%%% the standard forward difference quotient,
%%% u'(x_j) (approx)= [u(x_{j+1}) - u(x_j)]/h.
%%% Try it and see the improvement for yourself!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D2 = diag(ones(N-1,1),-1) -2*diag(ones(N,1)) + diag(ones(N-1,1),1);

x = (a:h:b)'; % the vector [x_0 ; x_1 ; ... ; x_(N+1)]
xshort = x(2:end-1); % the vector [x_1 ; x_2 ; ... ; x_N]
Xinv = diag(1./xshort);

L = D2/h^2 + Xinv*D1/h + eye(N);  %% discrete form of the differential operator
                               %%   d^2/dx^2 + (1/x)d/dx + 1

f = xshort.^2; % the vector [f(x1) ; f(x2) ; ... ; f(xN)]

u = L\f; % approximate solution


% uactual is the actual solution to the BVP (thanks, Maple!)
uactual = 12*(  besselj(0, xshort)*(8*bessely(0, 20)-33*bessely(0, 10)) ...
              - bessely(0, xshort)*(-33*besselj(0, 10)+8*besselj(0, 20)) ) ...
              / (besselj(0, 20)*bessely(0, 10)-bessely(0, 20)*besselj(0, 10)) ...
              - 4 + xshort.^2;

figure(1),clf
plot(xshort,[u,uactual])