function rmserror=nma_fem(a,b,c,f,x0,u0,beta,len,exact_sol,nNodes)
% Solves numerically a second order ODE with constant coeff.
% using the symmetric Garlekin Finite Elements Methods.
%
% see the file nma_fem_driver.m on how to call this function.
%
% Solve a u''(t) + b u'(t) + c u(t) = f(t)
%
% with t over the range t0 to len
% and with initial conditions u(0)=u0
% and with u'(len)=beta

%by Nasser Abbasi. Sept 26,2006.

    xc=linspace(x0,len,nNodes);

    % This plots the shape functions.
    %    x=linspace(x0,len,1000);
    %    y=linspace(x0,len,1000);
    %    for i=1:nNodes
    %        for j=1:length(y)
    %            [v,d]=phi(i,x(j),nNodes,x0,len,xc);
    %            y(j)=v;
    %        end
    %        plot(x,y);
    %        hold on;
    %    end


    A    = build_stiffness_matrix(nNodes,xc,a,b,c);
    load = build_load_vector(a,u0,nNodes,xc,f,beta);

    % Now remove the first row and column from the stiffness matrix
    
    A(2,1)  = A(2,1)*u0;
    load(2) = load(2)-A(2,1);
    A       = A(2:end,2:end);
    load    = load(2:end);
    
    % SOLVE for unknowns
    
    q = A\load;
    q = [u0;q];

    % Plot the solution
    y = zeros(length(xc),1);
    for i=1:length(xc)
        y(i)=trial(xc(i),x0,len,xc,nNodes,q);
    end

    plot(xc,y,'ro');
    hold on;
    line(xc,y);
    
    % Calculate RMSerror. Use 50 points. Should be more than enough.
    NPOINTS = 50;
    x = linspace(x0,len,NPOINTS);
    rmserror = 0;
    for i = 1:length(x)
        y = trial(x(i),x0,len,xc,nNodes,q);
        t=x(i);  % USED for subs below. do not remove.
        yexact   = real(double(subs(exact_sol)));
        rmserror = rmserror+(y-yexact)^2;
    end

    rmserror = rmserror/NPOINTS;
    rmserror = sqrt(rmserror);   
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v = trial(x,x0,len,xc,nNodes,q)
    if x<x0||x>len
        error('in Trial. x outside range');
    end
    v = 0;
    for i = 1:nNodes
        [s,d] = phi(i,x,nNodes,x0,len,xc); %notice ignore d here.
        v     = v+ q(i)*s;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [v,d]=phi(i,x,nNodes,x0,len,xc)
    if(i<1 || i>nNodes)
        error('node number outside range');
    end
    if(x<x0 || x>len)
        error('x outside range');
    end

    h = xc(2)-xc(1);

    if i==1 
        if(x>xc(2))
            v = 0;
            d = 0;
        else
            v = (xc(2)-x)/h;
            d = -1/h;
        end
        return;
    end

    if i==nNodes
        if(x<xc(nNodes-1))
            v = 0;
            d = 0;
        else
            v = (x-xc(nNodes-1))/h;
            d = 1/h;
        end
        return;
    end

    if(x>xc(i+1) || x<xc(i-1) )
        v = 0;
        d = 0;
    else
        if(x<=xc(i))
            v = (x-xc(i-1))/h;
            d = 1/h;
        else
            v = (xc(i+1)-x)/h;
            d = -1/h;
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SEE my report for description of this function.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A=build_stiffness_matrix(nNodes,xc,a,b,c)
    
    A = zeros(nNodes,nNodes);
    h = xc(2)-xc(1);
    for i = 1:nNodes
        if i==1
            A(1,1) = -a/h+ a*xc(1)/h^2 - b/h +c*h/3 - a*xc(2)/h^2;
            A(1,2) = -a*xc(1)/h^2 +b/2 -c*h/6 + a*xc(2)/h^2;
        else
            if i==nNodes
                A(nNodes,nNodes-1) = -a*xc(nNodes-1)/h^2 -b/2 +c*h/6 +a*xc(nNodes)/h^2;
                A(nNodes,nNodes)   = a*xc(nNodes-1)/h^2 +b/2 +c*h/3 -a*xc(nNodes)/h^2;
            else
                A(i,i-1) = -a*xc(i-1)/h^2 -b/2 + c*h/6 + a*xc(i)/h^2;
                A(i,i)   = a*xc(i-1)/h^2 + 2*c*h/3 -a*xc(i+1)/h^2;
                A(i,i+1) = -a*xc(i)/h^2 + b/2 + c*h/6 + a*xc(i+1)/h^2;
            end
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%Called to integrate the 'force' function
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v=integrand(x,f,xj)
    v = (x-xj).*f(x);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SEE my report for more description of this
% This build the load vector.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function load = build_load_vector(a,u0,nNodes,xc,f,beta)

    load = zeros(nNodes,1);
    h    = xc(2)-xc(1);
  
        
    load(1) = -(1/h)* quad(@(x)integrand(x,f,xc(2)),xc(1),xc(2)) ;
    load(nNodes) = (1/h)* quad(@(x)integrand(x,f,xc(nNodes-1)),xc(nNodes-1),xc(nNodes)) - a*beta;
    for i = 2:nNodes-1
        load(i) = quad(@(x)integrand(x,f,xc(i-1)),xc(i-1),xc(i)) - ...
                  quad(@(x)integrand(x,f,xc(i+1)),xc(i),xc(i+1));
        load(i) = (1/h)*load(i) ;
    end
end