function [ ErrL2, ErrH1 ] = solveFEM_assign1( N )

% solves the problem u_{xx} = - pi^2 * sin( x ) for x \in (-1,1) 
% u(-1) = u(1) = 0
% we use FEM, with polynomials
% phi_1 = (x-1)*(x+1) ... phi_i = x^(i-1) * (x-1)*(x+1) ...
% the input N is the number of polynomials that we can use
% the function returns the error in L2 and H1 norms
% integration is done symbolically


syms x; % this activates the symbolic library

% assign matrix T and vector R
T = zeros( N, N );
R = zeros( N, 1 );

% add the entries for T and R
for k = 1:N
    phi_k = x^(k-1)*(x-1)*(x+1); % get the k-th polynomial
    for i = 1:N
        phi_i = x^(i-1)*(x-1)*(x+1); % get the i-th polynomial
        % add the entries to the matrix
        T(k,i) = solveFEM_assign1_B( phi_i, phi_k );
    end;
    R(k) = solveFEM_assign1_R( phi_k ); % form the right-hand side
end;

% solve the projected finita dimensional problem
disp([' Condition number of the matrix is ',num2str(cond( T ))]);
C = T \ R;

% use the vecotr C to form back the solution, sum_{i=0}^N c_i phi_i
u = 0;
for i = 1:N
    phi_i = x^(i-1)*(x-1)*(x+1); % get the i-th test function
    u = u + C(i) * phi_i; % add C(i) times the i-th test function
end;

% Plot the solution
xx = -1:0.05:1; % pick points on the x-axis
pp = zeros( size( xx ) ); % this would be the computer approximation
uu = zeros( size( xx ) ); % this would be the exact solution
for i = 1:max(size(xx))
    pp(i) = subs( u, xx(i) ); % get the computer approximation at xx(i)
    uu(i) = subs( sin(pi*x), xx(i) ); % get the exact solution at xx(i)
end;
plot( xx, pp, xx, uu ); % plot

% get the errors
ErrL2 = solveFEM_assign1_L2( u );
ErrH1 = solveFEM_assign1_H1( u );

function [ I ] = solveFEM_assign1_B( u, phi )
% int( f, x, -1, 1 ) integrates a function f with respect to x from -1 to 1
% diff( f, x ) differentiates f with respect to x
% eval converts a symbolic value to a numerical value (i.e. pi becomes 3.14)
syms x;
I = - eval( int( diff( u, x ) * diff( phi, x ), x, -1, 1 ) );

function [ I ] = solveFEM_assign1_R( phi )
syms x;
I = eval( int( - pi^2 * sin( pi*x ) * phi, x, -1, 1 ) );

function [ E ] = solveFEM_assign1_L2( u )
syms x;
E = sqrt( eval( int( ( u - sin( pi * x ) )^2, x, -1, 1 ) ) );

function [ E ] = solveFEM_assign1_H1( u )
syms x;
E = sqrt( eval( int( ( u - sin( pi * x ) )^2 + ( diff( u, x ) - pi*cos( pi * x ) )^2, x, -1, 1 ) ) );

