%
% We solve the BVP
%     ep*y''+c1*y'+c2*y=c3 
%     y(a)=ya,  y(b)=yb
% The numerical solution is stored in two vectors x and y
%
clear
ep=0.01;  % the small parameter
%
a=0;	% one boundary condition is imposed at x=a
b=1;	% second boundary condition is imposed at x=b
ya=1; % boundary value at x=a
yb=-1;  % boundary value at x=b
N=2^11;
dx=(b-a)/N; % spatial step size
N1=N-1;
%
% We discretize the BVP as a large sparse linear system 
%
x=a+[0:N]'*dx;
c1=-(x(2:N)-1/3);
c2=-ones(N1,1);
c3=zeros(N1,1);;
%
ii=zeros(N1*3,1);
jj=zeros(N1*3,1);
ss=zeros(N1*3,1);
%
b=c3;
b(1)=b(1)-ya*ep/dx^2+ya*c1(1)/(2*dx);
b(N1)=b(N1)-yb*ep/dx^2-yb*c1(N1)/(2*dx);
ic=0;
%
ii(ic+1:ic+N1)=[1:N1]';
jj(ic+1:ic+N1)=[1:N1]';
ss(ic+1:ic+N1)=c2-ep*2/dx^2*ones(N1,1);
ic=ic+N1;
%
ii(ic+1:ic+N1-1)=[2:N1]';
jj(ic+1:ic+N1-1)=[1:N1-1]';
ss(ic+1:ic+N1-1)=ep/dx^2*ones(N1-1,1)-c1(2:N1)/(2*dx);
ic=ic+N1-1;
%
ii(ic+1:ic+N1-1)=[1:N1-1]';
jj(ic+1:ic+N1-1)=[2:N1]';
ss(ic+1:ic+N1-1)=ep/dx^2*ones(N1-1,1)+c1(1:N1-1)/(2*dx);
ic=ic+N1-1;
%
A=sparse(ii(1:ic),jj(1:ic),ss(1:ic),N1,N1);
%
yI=A\b;
y=[ya; yI; yb];
%
% Again, the numerical solution is stored in two vectors x and y
%
figure(1)
clf reset
% We plot one in every 8 points
x2=x(1:8:N+1);
y2=y(1:8:N+1);
plot(x2,y2,'b-*')
set(gca,'fontsize',16)
%
% save data_01 x2 y2
%
