%y''=p(x)y'+q(x)y+r(x),a<=x<=b
%y(a)=af,y(b)=bt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%實例%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  p=inline('-2/x');
  q=inline('2/x^2');
  r=inline('sin(log10(x)/log10(exp(1)))/x^2');
  N=9;
  a0=1;b0=2;
  af=1;bt=2;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  h=(b0-a0)/(N+1);
  x=a0+h;
  a(1)=2+h*h*q(x);
  b(1)=-1+(h/2)*p(x);
  d(1)=-h*h*r(x)+(1+(h/2)*p(x))*af;
    for i=2:N-1
        x=a0+i*h;
        a(i)=2+h*h*q(x);
        b(i)=-1+(h/2)*p(x);
        c(i)=-1-(h/2)*p(x);
        d(i)=-h*h*r(x);
     end
     x=b0-h;
     a(N)=2+h*h*q(x);
     c(N)=-1-(h/2)*p(x);
     d(N)=-h*h*r(x)+(1-(h/2)*p(x))*bt;
     %%%%%%%%%追趕法%%%%%%%%%%%%%%%%%%
     %y=trisys(c,a,b,d)
     L(1)=a(1);
     u(1)=b(1)/a(1);
     for i=2:N-1
         L(i)=a(i)-c(i)*u(i-1);
         u(i)=b(i)/L(i);
     end
     L(N)=a(N)-c(N)*u(N-1);
     z(1)=d(1)/L(1);
     for i=2:N
         z(i)=(d(i)-c(i)*z(i-1))/L(i);
     end
     y(N)=z(N);
     for i=N-1:-1:1
         y(i)=z(i)-u(i)*y(i+1);
     end
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Y=[af,y,bt];
 for i=1:N+2
      x=a0+(i-1)*h;
      zj(i)=1.1392070132*x-0.03920701320/x^2-3*sin(log10(x)/log10(exp(1)))/10-cos(log10(x)/log10(exp(1)))/10;
 end
disp('下面兩列分別是數值解和近似解');
 re=[Y'      zj']