function fd1d_singular (e,M,a,b,ua,ub)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%       FD1D
%%       for MTH 453/553
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This function shows how to set up 
%% the finite difference approximation of a model 
%% two-point boundary value problem -eu''+u=f
%% The user must code expression for utrue and for f (rhsfun)
%% and set up /change boundary conditions 
%% default discretization is 0....M, with M a parameter
%% but since MATLAB can't number vectors from 0, we use 1..M+1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% M. Peszynska, 01/2017
%% © Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
set(0,'defaultlinelinewidth',2);
%%%%%%%%%% setup the true solution. (YOU MUST FIND IT AND CODE IT)
utrue=@(x)(x);
%%%%%%%%%%%%%%%%%%%%%%%% 
function v = rhsfun(x)
    v = 1+0*x;
end
%%%%%%%%%%%%%%%%%%%%%%%%%% use uniform grid
clf;
xplot = linspace(a,b);
plot(xplot,utrue(xplot)); axis([0 1 -0.1 1.5]);pause;
%%%%%%%%%%%%%%%%%%%%%%%%% 
dx = (b-a)/M;           
x = (a:dx:b)';
%%%%%%%%%%%%%%%%%%%%%%%%%
%% set up the matrix and the right hand side vector
f = rhsfun(x);
%%
A = sparse(M+1,M+1);     
for j=2:M
    %% set up jth row of the matrix
    A(j,j) = 2;
    A(j,j-1) = -1; 
    A(j,j+1) = -1; 
end
A = A / (dx*dx);
%%
A = e*A + eye(M+1,M+1);

%% apply the Dirichlet boundary conditions to the matrix 
%% and to the right hand side
A (1,1) = 1;
A (M+1,M+1) = 1; 
f (1,1) = ua;
f (M+1,1) = ub;

%% solve the linear system;
u = A\ f;

%% plot both the analytical and the numerical solutions
plot(xplot,utrue(xplot),'k',x,utrue(x),'bo',x,u,'r*');
h=legend('utrue','u','u-FD');
set(h,'FontSize',12);
%% compute the error vector
error = utrue(x) - u;
%% report the error in max norm
fprintf('%g    %g %g\n',dx,norm(error,inf)); axis([0 1 -0.1 1.5]);
title(sprintf('True and numerical solution, eps=%g dx=%g',e,dx),'FontSize',12);
end
        
    
