%{
      Test: from the FD solution of y'' + kappa^2 y = 0 on (0,1)
        with boundary conditions y(0) = 0, y'(1) = 1.
       Second-order central FD for interior nodes;
       second-order, one-side FD (Gear) for right-most node.
%}

clear all; close all; clc;

bc  = [0 1];

N = 100;
h = 1/N;
x = 0:h:1;
Kappa = 2.5;
kappa = Kappa*h;

Neq = N-1;
d = (-2+kappa^2)*ones(Neq,1);
l = ones(Neq,1);
u = ones(Neq,1);
b = zeros(Neq,1);

% Modification due to Dirichlet boundary conditions:
b(1)   = b(1)   - l(1)*bc(1);

d(end) = d(end)+4/3; l(end) = l(end)-1/3; b(end) = b(end)-(2/3)*h*bc(2);


y = thomas(l,d,u,b);

y = [0;y;(2*h*bc(2)-y(end-1)+4*y(end))/3];

if Kappa==0
    yan = x;
else
    yan = sin(Kappa*x)./(Kappa*cos(Kappa));
end

figure
plot(x,yan,'k-');
hold on
plot(x,y,'rx');
box on