%本程序是為了求解變系數奇異問題的四次樣條函數方法，其精度能達到六階，是目前最好結果！

clear
format long
e=1/10^4;
h=1/512;
nh=1/h;
y(1)=0;
y(nh+1)=0;
for i=1:nh+1
        x=(i-1)*h;
 % f(i)=-40*(x*(x^2-1)-2*e);
 f(i)=1+x*(1-x)+(2*sqrt(e)-x*(1-x)^2)*exp(-x/sqrt(e))+(2*sqrt(e)-x^2*(1-x))*exp(-(1-x)/sqrt(e));
  %p(i)=1+x;
  p(i)=1+x*(1-x);
    %    f(i)=1+2*sqrt(e)*(exp(-x/sqrt(e))+exp((x-1)/sqrt(e)))
   %f(i)=(i-1)*h;
end
for j = 1:nh-3
    if  j ~= nh-3
        A(j,j) =11*p(j+2)+51/2*e/h^2; 
        A(j+1,j) =2*p(j+2)-12*e/h^2;
        A(j+2,j)=-3*e/(4*h^2);
        A(j,j+1) =2*p(j+3)-12*e/h^2;
        A(j,j+2)=-3*e/(4*h^2);
    else
        A(j,j) =11*p(j+2)+51/2*e/h^2;
    end
end
A2=A(1:nh-3,1:nh-3);
A1(2:nh-2,2:nh-2)=A2;
A1(1,1)=25/2*p(2)+30*e/h^2;
A1(1,2)=5/4*p(3)-15*e/h^2;
A1(2,1)=2*p(2)-12*e/h^2;
A1(3,1)=-3*e/(4*h^2);
A1(nh-2,nh-1)=2*p(nh)-12*e/h^2;
A1(nh-3,nh-1)=-3*e/(4*h^2);
A1(nh-1,nh-2)=5/4*p(nh-1)-15*e/h^2;
A1(nh-1,nh-1)=25/2*p(nh)+30*e/h^2;
n=nh-3;
SM1=sparse(1:n,1:n,11*ones(1,n),n,n,n);
SM2=sparse(2:n,1:n-1,2*ones(1,n-1),n,n,n-1);
SM3=sparse(3:n,1:n-2,0*ones(1,n-2),n,n,n-2);
S1=SM1+SM2+SM2'+SM3+SM3';
B=full(S1);
    D(2:nh-2)=B*f(3:nh-1)';
    D(1)=5/4*(f(3)+10*f(2)+f(1))-y(1)*(5/4*p(1)-15*e/h^2);
    D(2)=D(2)+y(1)*(3*e/(4*h^2))+2*f(2);
    D(nh-2)=D(nh-2)+y(nh+1)*(3*e/(4*h^2))+2*f(nh);
    D(nh-1)=5/4*(f(nh-1)+10*f(nh)+f(nh+1))-y(nh+1)*(5/4*p(nh+1)-15*e/h^2);
    y(2:nh)= lufact(A1,D');
    for i=1:nh+1
            x=(i-1)*h;
   % u(i)=(exp(-(1-(i-1)*h)/e^(1/2))+exp(-(i-1)*h/e^(1/2)))/(1+exp(-1/e^(1/2)))-cos(pi*(i-1)*h)*cos(pi*(i-1)*h);
   % u(i)=1-(1-x)*exp(-x/sqrt(e))-x*exp((x-1)/sqrt(e));
      %u(i)=exp(-(i-1)*h/e^(1/2))+(i-1)*h;
%u(i)=40*x*(1-x);
u(i)=1+(x-1)*exp(-x/sqrt(e))-x*exp(-(1-x)/sqrt(e));
end
g1=y-u;
max(abs(g1))
%x=0:h:1;
%plot(x,y,'+',x,u,'-.')
%q=[x',y',u']