% Solve a Two-point ODE Boundary Value Problem (BVP) with Finite Difference Method (FDM)
% u''(x)-4u(x)=0,\quad u(0)=1, u(1)=3; (see Example 7.8 in page 277)
% Using ODE theory (method of undetermined coefficients), we can find its exact solution:
% u(x)=c_1e^{2x}+c_2e^{-2x} with c_1=\frac{3-e^{-2}}{e^2-e^{-2}}, c_2=\frac{e^2-3}{e^2-e^{-2}}.

c1=(3-exp(-2))/(exp(2)-exp(-2));c2=(exp(2)-3)/(exp(2)-exp(-2));
ufun=@(x) c1*exp(2*x)+c2*(exp(-2*x));%exact solution
ffun=@(x) zeros(size(x)); %right-hand size

% Following our Lecture of Section 7.2, we first define a uniform mesh of [0,1]: x_j=jh,  j=0,1,\cdots,m+1, with\ h=\frac{1}{m+1}.
% Here we divide it into (m+1) sub-intervals.

m=5; h=1/(m+1);%step size
x=linspace(0,1,m+2)'; %divide [0,1] into (m+2) points
figure(1)
plot(x,ufun(x),'o--'); %plot the exact solution
xlabel('x'); ylabel('u(x)');

% Then we will discretize the ODE BVP by the Central FDM to get a sparse linear system.
% Since the ODE itself holds at each interior point x_j, j=1,\cdots,m,which gives
% u''(x_j)-4u(x_j)=0, j=1,2,\cdots,m, with the boundary conditions u(x_0)=u(0)=1; u(x_{m+1})=u(1)=3.
% Now, using the central finite difference formula for approximating u''(x_j) to get
% u''(x_j)=\frac{u(x_{j-1})-2u(x_j)+u(x_{j+1})}{h^2}+O(h^2).
% Let u_j=u(x_j) and f_j=f(x_j), then the above ODE at each point now reads
% \frac{u_{j-1}-2u_j+u_{j+1}}{h^2}+{\bf O(h^2)}-4u_j=0.
% Upon dropping the O(h^2) truncation error term, we get a linear system of finite difference formulas (see lecture for the matrix formulation pattern)
% \frac{U_{j-1}-2U_j+U_{j+1}}{h^2}-4U_j=0, j=1,2,\cdots,m,
% where  U_j\approx u_j denotes the finite difference numerical approximation.
% By formualting this into a sparse system of linear equations A_h \vec{U}_h=\vec{f}_h, where
%% A_h=\left[\matrix{ -2/h^2-4 & 1/h^2 & 0 & \cdots&\cdots & 0 \cr
% 1/h^2 &-2/h^2-4 &1/h^2 &0 &\cdots &0 \cr
% 0&1/h^2 &-2/h^2-4 &1/h^2 &\cdots &0\cr
% \vdots &\ddots &\ddots &\ddots &\ddots &\vdots \cr
% 0&\cdots&0&1/h^2 &-2/h^2-4 &1/h^2 \cr
% 0&0&\cdots&0&1/h^2 &-2/h^2-4}\right]
% \vec{f}_h=[0-1/h^2,0,\cdots,0,0-3/h^2]^T,
% we can solve the linear system for \vec{U}_h=[U_1,U_2,\cdots,U_m]^T, which gives the numerical solution.
% Under certain assuptions, we can show it converges to the exact solution in the sense that
% \| \vec{U}_h-\vec{u}\|_\infty=O(h^2)\to 0 as h\to 0, (other vector norms can also be used)
% where \vec{u}=[u_1,u_2,\cdots,u_m]^T is the exact solution at the same grid points.
% Below is the vectorized codes for building this system using sparse matrices (based on spdiags):

e=ones(m,1);
Ah=spdiags([e/h^2 -2*e/h^2-4 e/h^2],(-1:1),m,m);%tridiagonal matrix
fh=ffun(x(2:end-1));%define rhs over interior nodes
fh(1)=fh(1)-1/h^2;fh(end)=fh(end)-3/h^2;%enforce boundary condition

% Once the system is constructed, we can for now solve it using backslash solver (sparse direct solver)

Uh=Ah\fh;%Solve the system by backslash solver, finite different solution
uh=ufun(x(2:end-1));%exact solution over interior nodes
maxerr=norm(Uh-uh,"inf") %this error should of O(h^2)

% The computed maxerr is of O(h^2), which means it will be four time smaller if h is halved.
% Plot the finite different solution versus the exact solution, see almost not difference

figure(2);hold on; plot(x(2:end-1),uh,'o--'); %plot the exact solution
plot(x(2:end-1),Uh,'x--r') %plot the finite difference solution
title('Exact solution vs finite different solution')

% It is sometimes better to see their difference by plotting the point-wise errors (no errors at the boundary)

figure(3); plot(x(2:end-1),Uh-uh,'-s')
xlabel('x');title('point-wise error at interior nodes')

% Often, one can use loglog plot of maxerr as a function of h to verify the order of accuracy.
% Notice maxerr=C*h^2 implies log(maxerr)=log(C)+2*log(h), which is a log-log line with slope 2.
% Below is a summary of the above codes within a loop of refined mesh sizes to generate an error log-log plot,
% where you need to understand that h limits to zero when you look from right to the left.

figure(4)
mlist=2.^(5:10);
hlist=1./(mlist+1);
maxerr=zeros(size(mlist));
for k=1:length(mlist)
    m=mlist(k); h=1/(m+1); x=linspace(0,1,m+2)'; e=ones(m,1);
    Ah=spdiags([e/h^2 -2*e/h^2-4 e/h^2],(-1:1),m,m);%tridiagonal matrix
    fh=ffun(x(2:end-1));%define rhs over interior nodes
    fh(1)=fh(1)-1/h^2;fh(end)=fh(end)-3/h^2;%enforce boundary condition
    Uh=Ah\fh;%Solve the system by backslash solver, finite different solution
    uh=ufun(x(2:end-1));%exact solution over interior nodes
    maxerr(k)=norm(Uh-uh,"inf"); %this error should of O(h^2)
end
p = polyfit(log(hlist),log(maxerr),1); slope=p(1) %Fix by a line to check the slop
loglog(hlist,maxerr,'d-'); xlabel('h');ylabel('maxerr');%should give a straight line with slop 2 