function [ xi, u ] = quadfem ( load, N );
% QuadFEM
%
% This m-file is part of the material of the lecture
%   Scientific Computing and Software Tools
% at the Math Department of Halle University, Germany.
% Summer term 2004
%
% Author :      Prof. Dr. M. Arnold, arnold@mathematik.uni-halle.de
% Version of :  Jun 22, 2004
%
% Solution of boundary value problem by piecewise quadratic finite elements
%
% Parameters
%   load   (input)  : function handle to define the load f(x)
%                     Example: see myload.m
%   N      (input)  : number of inner grid points
%   xi     (output) : Lagrange interpolation points
%   u      (output) : FE solution vector
%
% Example
%   [ xi, u ] = quadfem ( @myload, 999 );
%   plot ( xi, u );

% -> define the grid
h  = 1/(N+1);
x  = linspace ( 0, 1, 1+N+1 );
xi = zeros ( 1+2*N+2, 1 );
xi(1:2:end) = x;
xi(2:2:end) = ( x(1:end-1) + x(2:end) ) / 2;

% -> function values of form functions at Gauss nodes
c = 0.5 + sqrt(3)/6 * [ 1 -1 ]';

phil = ( 1 - 2*c ) .* ( 1 - c );
phic = 4 * ( 1 - c ) .* c;
phir = ( 2*c - 1 ) .* c;

phipl = (  4*c - 3 ) / h;
phipc = ( -8*c + 4 ) / h;
phipr = (  4*c - 1 ) / h;

% -> get matrix elements by Gauss quadrature, s = 2
aoc = h * sum ( phipc .* phipc ) / 2;
aop = h * sum ( phipc .* phipl ) / 2;

aec  = h * sum ( phipl .* phipl + phipr .* phipr ) / 2;
aep  = aop;
aepp = h * sum ( phipl .* phipr ) / 2;

% -> FE matrix
n = 2*N + 1;

diagc = zeros ( n, 1 );
diagc(1:2:end) = aoc;
diagc(2:2:end) = aec;

diagp = zeros ( n, 1 );
diagp(1:2:end) = aep;
diagp(2:2:end) = aop;

diagpp = zeros ( n, 1 );
diagpp(2:2:end) = aepp;

a = spdiags ( [ diagpp, diagp ], -2:-1, n, n )  + ...
    spdiags ( [ diagpp, diagp ], -2:-1, n, n )' + ...
    spdiags ( diagc, 0, n, n );

% -> FE right hand side
b = zeros ( n, 1 );
f = zeros ( size ( c ) );

for l=0:N,
    i = 2*l + 1;
    for ic=1:length(c),
        xact = x(1+l) + c(ic)*h;
        f(ic) = load ( xact );
    end;
    b(i) = h * sum ( phic .* f ) / 2;
    if l>0,  b(i-1) = b(i-1) + h * sum ( phil .* f ) / 2;  end;
    if l<N,  b(i+1) = h * sum ( phir .* f ) / 2;  end;
end;

% -> FE solution
u = [ 0; a\b; 0 ];