function gunz()

% Linear finite element for sturm_liouville problem
% Originally by M. Gunzburger 4/1977
% Translated to matlab by greg
%
% -D(P(x)Du(x))+Q(x)u(x)=F(x), x in (a,b)
%         Du(a)=Du(b)=0
clear global xn;
clear global node;

he=zeros(100,1); xq=zeros(100,2); g=zeros(101,1); a=zeros(101,2); indx=zeros(101,1);
global xn node
xn=zeros(101,1); 
node=zeros(101,2);

%....This is a sample problem for a problem of the above type with
%....a=PI,b=3*PI, uniform mesh, with p, q, f, exact given below

%....Various funstions

grid = @(x) x;
p= @(x) 1;
q= @(x) 1;
f= @(x) 1.25*sin(0.5*x);
exact= @(x) sin(0.5*x);

%...Input

disp(sprintf('Choose one of the following decomposition for the system'));
disp(sprintf(' 1: Gaussian elimination (original algorithm)'));
disp(sprintf(' 2: Cholesky decomposition & backsolving'));
disp(sprintf(' 3: Tridiagonal LU decomposition algorithm & bachsolving'));
iflag=input('(Just type the number of the special algorithm)   ');

n=-1;
while n<=0
n = input('Give the number of intervals   ');
end
n=fix(n);


%....Endpoints of the interval (a,b)

% pi= 4*atan(1);
ax=pi;
bx=3*pi;

%....Construct grid (uniform mesh) and order unknows
h = (bx-ax)/n;
xx = ax-h;
np = n+1;
i = 0;
for in=1:np
    xx=xx+h;
    xn(in)=grid(xx);
    i=i+1;
    indx(in)=i;
end
nu = i;

%....Element information

for ie=1:n
    xl=xn(ie);
    in=ie+1;
    xr=xn(in);
    he(ie)=xr-xl;
    node(ie,1)=ie;
    node(ie,2)=in;
    xq(ie,1)=xl;
    xq(ie,2)=xr;
end

%....Matrix assembly by elements

for i=1:nu
    g(i)=0;
    a(i,1)=0;
    a(i,2)=0;
end
for ie=1:n
    h=he(ie)*.5;
    for iq=1:2
        x=xq(ie, iq);
        for ic=1:2
            ii= node(ie, ic);
            bb= rofu(x, ie, ic, 1);
            bx= rofu(x, ie, ic, 2);
            i= indx(ii);
            g(i)=g(i)+bb*f(x)*h;
            for jc=1:2
                jj= node(ie, jc);
                bbb= rofu(x, ie, jc, 1);
                bbx= rofu(x, ie, jc, 2);
                j=indx(jj);
                if (j > i)
                    continue
                end
                aij=h*(p(x)*bx*bbx + q(x)*bb*bbb);
                ij=2-i+j;
                a(i,ij)=aij+a(i,ij);
            end
        end
    end
end

%....Compute the relative machine precision
%....(machine epsilon, see forsythe, malcon & moler)

% eps=1;
% while 1+eps>1
%     eps=0.5*eps;
% end

%USING eps INSTEAD

%...Solve algebraic system AQ=G

switch iflag
    case 1
        %..Gaussian elimination
        for i=2:nu
            if(abs(a(i-1,2))<=eps)
                error('Gauss elimination not possible');
            end
            a(i,1)=a(i,1)/a(i-1,2);
            a(i,2)=a(i,2)-a(i-1,2)*a(i,1)*a(i,1);
            g(i)=g(i) - g(i-1)*a(i,1);
        end
        if(abs(a(nu,2))<=eps)
            error('Gauss elimination not possible');
        end
        g(nu)=g(nu)/a(nu,2);
        for i=2:nu
            ic=nu-i+1;
            g(ic)=g(ic)/a(ic,2)-g(ic+1)*a(ic+1,1);
        end
         
    case 2
        %..Cholesky decomposition & backsolving
        eps2=eps*eps;
        if(a(1,2) > eps2)
            a(1,2) = sqrt(a(1,2));
            g(1)=g(1)/a(1,2);
        else
            error('Cholesky decomposition not possible');
        end
        for i=2:nu
            a(i,1)=a(i,1)/a(i-1,2);
            ww=a(i,2)-a(i,1)*a(i,1);
            if (ww <= eps2)
                error('Cholesky decomposition not possible');
            end
            a(i,2) = sqrt(ww);
            g(i)=(g(i)-g(i-1)*a(i,1))/a(i,2);
        end
        g(nu)=g(nu)/a(nu,2);
        for i=2:nu
            ic=nu-i+1;
            g(ic)=(g(ic)-g(ic+1)*a(ic+1,1))/a(ic,2);
        end
        
        
    case 3
        %..Tridiagonal L_U decomposition algorithm & backslving
        g(1)=g(1)/a(1,2);
        for i=2:nu
            xx=a(i,1)/a(i-1,2);
            a(i,2)=a(i,2)-a(i,1)*xx;
            if(abs(a(i,2)) <= eps)
                error('Tridiagonal L_G decomposition not posseble');
            end
            g(i)=(g(i)-g(i-1)*a(i,1))/a(i,2);
            a(i,1)=xx;
        end
        for i=2:nu
            ic=nu-i+1;
            g(ic)=g(ic)-g(ic+1)*a(ic+1,1);
        end
        
    otherwise
        error('error solving system AQ=G. Method unknow...');
end

%....Output solution

disp(sprintf('\nnumber of intervals=%d\nx\t\t\tsolution\terror\t\texact', n));

errmax=0;
errl2=0;
for i=1:nu
    x=xn(i);
    s=g(i);
    err=abs(s-exact(x));
    disp(sprintf('%f\t%f\t%f\t%f', x, s, err, exact(x)));
    if(i ~= n+1)
        errl2= errl2+he(i)*(exact(x)-s)*(exact(x)-s);
    end
    if (err <= errmax)
    continue
    end
    errmax=err;
end
        
errl2= sqrt(errl2);
disp(sprintf('\n\tMaximum error=\t%g\n\tDiscrete L2_error=\t%f', errmax, errl2));

%GOTO BEGIN

end




function rofu=rofu(x,ie, ic, j)
global xn node
ic1 = ic;
ic2 = mod(ic,2)+1;
i1= node(ie, ic1);
i2= node(ie, ic2);
x1= xn(i1);
x2= xn(i2);
if(j==2)
    rofu=1.0/(x1-x2);
else
    rofu=1.0*(x-x2)/(x1-x2);
end
end








