function gunz_cbspl()
% ***********************************************************************
%
%      CUBIC B_SPLINE FINITE ELEMENT PROGRAM FOR 2_POINT B.V.P.
%      OF THE FORM :
%
%         -(P*U')' + Q*U = F , X IN (A,B)
%
%           U'(A) = U(B) = 0   ( ' DENOTES DIFFERENTIATION W.R. TO X )
% ***********************************************************************

% -----VARIOUS FUNCTIONS
HPI = 1.5707963267948966D0;
P = @(x) 1.0D0;
Q = @(x) 1.0D0;
F = @(x) (1.0D0+HPI*HPI)*cos(HPI*x);
EXACT = @(x) cos(HPI*x);
EXACTX = @(x) -HPI*sin(HPI*x);

% -----INPUT
N = input('Give the number of intervals: ');
if N <= 0; error 'Not enough intervals'; end

% -----MEMORY ALLOCATION
A = zeros(N+3,4); G = zeros(N+3,1); XN = zeros(N+3,1); INDX = zeros(N+3,1);
NODE = zeros(N,4); XQ = zeros(N,3);
NDIM = N+3;
% % COMMON /A1/N /A2/HH

% -----ENDPOINTS OF THE INTERVAL
AX = 0.D0;
BX = 1.D0;

% -----WEIGHTS OF GAUSS-LEGENDRE INTEGRATION FORMULA ( 3 POINTS )
WQ(1) = (BX-AX)*5.0D0/1.8D01;
WQ(2) = (BX-AX)*8.0D0/1.8D01;
WQ(3) = WQ(1);
CON = sqrt(0.6D0);

% -----CONSTRUCT GRID AND ORDER UNKNOWNS
HH = N;
H = (BX-AX)/HH;
X = AX-2.0D0*H;
NP = N+3;
I = 0;
for IN=1:NP
    X = X+H;
    XN(IN) = X;
    I = I+1;
    INDX(IN) = I;
end
NU = NP-1;
INDX(NP) = 0;

% -----ELEMENT INFORMATION
for IE=1:N
    XL= XN(IE+1);
    XR= XN(IE+2);
    NODE(IE,1)= IE;
    NODE(IE,2)= IE+1;
    NODE(IE,3)= IE+2;
    NODE(IE,4)= IE+3;
    XQ(IE,2)= XL+XR;
    XQ(IE,1)= (XQ(IE,2)-H*CON)*0.5D0;
    XQ(IE,3)= XQ(IE,2)-XQ(IE,1);
    XQ(IE,2)= XQ(IE,2)*0.5D0;
end

% -----MATRIX ASSEMBLY BY ELEMENTS
for I=1:NU
    G(I) = 0.0D0;
    for J=1:4
        A(I,J) = 0.0D0;
    end
end
G(NP) = 0.0D0;
for IE=1:N
    for IQ=1:3
        X = XQ(IE,IQ);
        for IC=1:4
            II = NODE(IE,IC);
            I = INDX(II);
            if( I == 0 ); continue; end
            BB = BSPL(X,IE,IC,1, N,HH);
            BX = BSPL(X,IE,IC,2, N,HH);
            G(I) = G(I) + BB*F(X)*WQ(IQ)*H;
            for JC=1:4
                JJ = NODE(IE,JC);
                J = INDX(JJ);
                if( J > I ); continue; end
                if( J == 0 ); continue; end
                BBB = BSPL(X,IE,JC,1, N,HH);
                BBX = BSPL(X,IE,JC,2, N,HH);
                AIJ = (P(X)*BX*BBX + Q(X)*BB*BBB) * H * WQ(IQ);
                IJ = 4 - I + J;
                A(I,IJ) = A(I,IJ) + AIJ;
            end
        end
    end
end

% -----SOLVE ALGEBRAIC SYSTEM
A = FACTOR(NDIM,NU,3,A);
G = SOLVE(NDIM,NU,3,A,G);
 
% -----OUTPUT SOLUTION
NP1= NP - 1;
%   70  FORMAT(/1X,'NUMBER OF INTERVALS = ',I3//)
fprintf('Number of intervals = %d\n', N);
%   80  FORMAT(8X,' X',18X,'EXACT',15X,'APPROXIMATE'/)
fprintf('X\t\t\t\tEXACT\t\t\tAPPROXIMATE\n');
for I=2:NP1
    X= XN(I);
    TEMP= EXACT(X);
    TEMP1= 0.25D0 * ( G(I-1) + G(I+1) ) + G(I);
    if( I == NP1 ); TEMP1= 0.0D0; end
    % 100  FORMAT(D13.4,2D21.8)
    fprintf('%e\t%e\t%e\n', X, TEMP, TEMP1); %warning
end

% -----COMPUTE ERRORS
EL2= 0.0D0;
EH1= 0.0D0;
for IE=1:N
    for IQ=1:3
        X= XQ(IE,IQ);
        SOL= 0.0D0;
        SOLX= 0.0D0;
        for IC=1:4
            II= NODE(IE,IC);
            I= INDX(II);
            if( I == 0 ); continue; end
            BB= BSPL(X,IE,IC,1, N,HH);
            BX= BSPL(X,IE,IC,2, N,HH);
            TEMP= G( IE + IC - 1 );
            SOL= SOL + BB * TEMP;
            SOLX= SOLX + BX * TEMP;
        end
        TEMP= SOL - EXACT(X);
        TEMP1= SOLX - EXACTX(X);
        EL2= EL2 + TEMP * TEMP * WQ(IQ) * H;
        EH1= EH1 + TEMP1 * TEMP1 * WQ(IQ) * H;
    end
end
EH1= sqrt( EH1 + EL2 );
EL2= sqrt( EL2 );
%  130  FORMAT(//8X,' H',16X,'L_2 ERROR',15X,'H-1 ERROR'/)
fprintf('\nH\t\t\t\tL2 ERROR\t\tH1 ERROR\n');
%  140  FORMAT(D15.4,2D21.8)
fprintf('%e\t%e\t%e\n', H, EL2, EH1);

end


%----------------BSPL----------------
function BSPL = BSPL(X,IE,I,J, N, HH)

Y= X*HH - (IE-1);
if( ~(IE==N && I>1) ) %---b.c.---
    if( J==1 ) %---deriv---
        switch(I)
            case 1
                CON = 1.0D0 - Y;
                BSPL = CON*CON*CON*0.25D0;
                return
            case 2
                CON = 1.0D0 - Y;
                BSPL = 0.25D0*( 1.0D0 + 3.0D0*CON*( 1.0D0 + CON*(1.0D0-CON) ) );
                return
            case 3
                BSPL = 0.25D0*( 1.0D0 + 3.0D0*Y*( 1.0D0 + Y*(1.0D0-Y) ) );
                return
            case 4
                BSPL = 0.25D0*Y*Y*Y;
                return
        end
    else %---deriv---
        switch(I)
            case 1
                CON = 1.0D0 - Y;
                BSPL = -0.75D0*CON*CON*HH;
                return
            case 2
                CON = 1.0D0 - Y;
                BSPL = -0.75D0*( 1.0D0 + CON*( 2.0D0 - 3.0D0*CON ) )*HH;
                return
            case 3
                BSPL = 0.75D0*( 1.0D0 + Y*( 2.0D0 - 3.0D0*Y ) )*HH;
                return
            case 4
                BSPL = 0.75D0*Y*Y*HH;
                return
        end
    end %---deriv---
    
else %---b.c.---
    
    if( J == 1 )
        CON = 0.25D0*Y*Y*Y;
        CON1 = 1.0D0 - Y;
        if( I == 2 ); BSPL = 0.25D0*( 1.0D0 + 3.0D0*CON1*( 1.0D0 +  CON1*(1.0D0-CON1) ) ) - CON;  end
        if( I == 3 ); BSPL = 0.25D0*( 1.0D0 + 3.0D0*Y*( 1.0D0 +  Y*(1.0D0-Y) ) ) - 4.0D0*CON; end
        return
    else
        CON = 0.75D0*Y*Y*HH;
        CON1 = 1.0D0 - Y;
        if( I == 2 ); BSPL = -0.75D0*( 1.0D0 + CON1*( 2.0D0 -3.0D0*CON1 ) )*HH - CON; end
        if( I == 3 ); BSPL = 0.75D0*( 1.0D0 + Y*( 2.0D0 -3.0D0*Y ) )*HH - 4.0D0*CON; end
        return
    end
end %---b.c.---

end

% ----------------FACTOR----------------
function A = FACTOR(NDIM,N,IB,A)

IB1 = IB + 1;
for I=1:N
    SUM = 0.0D0;
    if(I~=1)
        IM1 = I - 1;
        K0 = max([1,IB1-I+1]);
        II = I;
        JJ = IB1;
        for K=K0:IB
            II = II-1;
            JJ = JJ-1;
            SUM = SUM + A(II,IB1)*A(I,JJ)*A(I,JJ);
        end
    end
    A(I,IB1) = A(I,IB1) - SUM;
    PIV = 1.0D0/A(I,IB1);
    if (I==N); return; end
    IP1 = I+1;
    J0 = min([N,I+IB]);
    for J=IP1:J0
        JI = IB1+I-J;
        SUM = 0.0D0;
        if (I==1); A(J,JI) = (A(J,JI)-SUM)*PIV; 
            continue; end
        II = I;
        JJ = JI;
        IJ = IB1;
        K0 = max([1,I-IB,J-IB]);
        if (K0>IM1); A(J,JI) = (A(J,JI)-SUM)*PIV; 
            continue; end
        for K=K0:IM1
            II = II-1;
            IJ = IJ-1;
            JJ = JJ-1;
            SUM = SUM + A(II,IB1)*A(I,IJ)*A(J,JJ);
        end
        A(J,JI) = (A(J,JI)-SUM)*PIV;
    end
end
end


% ----------------SOLVE----------------
function G = SOLVE(NDIM,N,IB,A,G)

IB1 = IB + 1;
G(1)= G(1)/A(1,IB1);
for I=2:N
    SUM = 0.0D0;
    K0 = max([1,IB1-I+1]);
    for J=K0:IB
        II = I - IB1 + J;
        SUM = SUM + A(I,J)*G(II)*A(II,IB1);
    end
    G(I) = (G(I)-SUM)/A(I,IB1);
end
NM1 = N - 1;
for K=1:NM1
    I = N - K;
    IP1 = I + 1;
    SUM = 0.0D0;
    K0 = min(N,I+IB);
    JJ = IB1;
    II = I;
    for J=IP1:K0
        II = II + 1;
        JJ = JJ - 1;
        SUM = SUM + G(J)*A(II,JJ);
    end
    G(I) = G(I) - SUM;
end
end






