%{
      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);
b(end) = b(end) - u(end)*bc(2);

%% d(N-1) = 2 + kappa^2; l(N-1) = 0; b(N-1) = -2*h;
%%d(N-1) = -1 + kappa^2; l(N-1) = 1; b(N-1) = -h;

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

%%A = spdiags([[l(2:end);nan],d,u],[-1,0,1],N-1,N-1);
% % % A = zeros(Neq,Neq);
% % % for j=1:Neq
% % %   A(j,j) = d(j);
% % %   if (j>1) A(j,j-1) = l(j); end
% % %   if (j<Neq) A(j,j+1) = u(j); end  
% % % end
% % % 
% % % y = A\b;

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

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

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