%求解一維問題的有限元法的通用子程序
%輸入變量解釋：a,b為變量的初末值，n為等分的份數，mo,mn為函數的初末值
function u=fem(a,b,n,m0,mn)
h=(b-a)/n;
%用MATLAB的內部函數來求解有限元方程系數的積分，b的求解是用外部定義的子函數求解
syms h1 kesi1;
f1=1/h+h*(1-kesi1)^2;
f1=int(f1,kesi1,0,1);
ai1i1=eval(f1);
syms h2 kesi2;
f2=-1/h+h*kesi2*(1-kesi2);
f2=int(f2,kesi2,0,1);
h2=h;
ai1i=eval(f2);
syms h3 kesi3;
f3=1/h+h*kesi3^2;
f3=int(f3,kesi3,0,1);
h3=h;
aii=eval(f3);
%求解有限元方程的系數矩陣和右端的常數向量
A=zeros(n+1,n+1);
B=zeros(n+1);
for i=2:n+1
    A1=zeros(n+1,n+1);
    B1=zeros(n+1);
    A1(i-1,i-1)=ai1i1;
    A1(i-1,i)=ai1i;
    A1(i,i-1)=ai1i;
    A1(i,i)=aii;
    A=A+A1;
    x(i-1)=h*(i-2);
    B1(i-1)=b1(x(i-1),h);
    B1(i)=b2(x(i-1),h);
    B=B+B1;
end
%調整有限元方程的系數矩陣和右端的常數向量
A=A(2:n,2:n);
B=B(2:n);
x=x(2:n);
%用不選主元的gauss消去法求解有限元方程
[L,U]=nongauss(A);
y=downsan(L,B');
u=upsan(U,y);
%把函數值變為變換前的函數值
u=u+x'*(exp(1)-exp(-1));