
% parameter of ODE 
b = @(x) ones(size(x,1),1);
%b = @(x) -(x-0.5);
%b = @(x) cos(2*pi*x);
c = @(x) zeros(size(x,1),1);
f = @(x) ones(size(x,1),1);

epsilon = 0.01;

% boundary conditions 
alpha0 = 1;
beta0 = 0;
gamma0 = 0;
alpha1 = 1;
beta1 = 0;
gamma1 = 0;


% parameter for discretization
hunif=0.005;
mesh = (0:hunif:1)';

% non uniform mesh
%mesh = [(0:0.1:0.9)';(0.91:0.01:1)'];

% compute stiffness matrix
N=size(mesh,1)-1;
h = mesh(2:end)-mesh(1:(end-1));
hj = h(2:end);
hjm1 = h(1:(end-1));

Ajj = -2./(hjm1.*(hjm1+hj));
Ajjp1 = 2./(hj.*hjm1);
Ajjp2 = -2./(hj.*(hj+hjm1));
A = spdiags([Ajj,Ajjp1,Ajjp2],[0,1,2],N-1,N+1);

Bjjp1 = -b(mesh(2:(end-1)))./hj;
Bjjp2 = -Bjjp1;
Bplus = spdiags([Bjjp1,Bjjp2],[1,2],N-1,N+1);

Bjj = -b(mesh(2:(end-1)))./hjm1;
Bjjp1 = -Bjj;
Bminus = spdiags([Bjj,Bjjp1],[0,1],N-1,N+1);

Bjj = -(hj./hjm1).*b(mesh(2:(end-1)))./(hj+hjm1);
Bjjp1 = -((hjm1./hj)-(hj./hjm1)).*b(mesh(2:(end-1)))./(hj+hjm1);
Bjjp2 = (hjm1./hj).*b(mesh(2:(end-1)))./(hjm1+hj);
B0 = spdiags([Bjj,Bjjp1,Bjjp2],[0,1,2],N-1,N+1);

C = spdiags(c(mesh(2:(end-1))),[1],N-1,N+1);

% change for different method
B = B0;

S = epsilon*A+B+C;

F = f(mesh(2:(end-1)));


% include boundary conditions
% alpha0 u(0) + beta0 u'(0) = gamma0
% alpha1 u(1) + beta1 u'(1) = gamma1

Sbc = [alpha0 - beta0/h(1), beta0/h(1), sparse(1,N-1);
       S;
       sparse(1,N-1), -beta1/h(end), alpha1 + beta1/h(end)];

Fbc = [gamma0;F;gamma1];

% solve linear system
u = Sbc\Fbc;

figure; 
plot(mesh,u);

% % plot first derivative
% du = (u(2:end)-u(1:(end-1)))./h;
% figure;
% plot((mesh(1:(end-1))+mesh(2:end))/2,du);
% 
% % plot second derivative
% d2u = -A*u;
% figure;
% plot(mesh(2:(end-1)),d2u);


