function solveFEM_assign3_nonuniform( N, break_point )

% get the mesh
%[ mesh_points, mesh_index ] = solveFEM_assign3_formUniformMesh( N );
[ mesh_points, mesh_index ] = solveFEM_assign3_formNonUniformMesh( N, break_point );

%solveFEM_assign3_B( @solveFEM_assign3_phi_x_L, @solveFEM_assign3_phi_x_L, @solveFEM_assign3_phi_x_L, 0, 1 )

T = sparse( N-1, N-1 );

% cycle over the elements to compute all the integrals for the left hand side
for e = 1:N
    
    [ w, x, xi ] = solveFEM_assign3_gauss3point( mesh_points(e,1), mesh_points(e,2) );
    
    iL = mesh_index(e,1); % get the left and right indexes
    iR = mesh_index(e,2);
    
    inv_dx = 2 / ( mesh_points(e,2) - mesh_points(e,1) );
    
    % check the left test function
    if ( iL > 0 )
        % left to left
        T( iL, iL ) = T( iL, iL ) + w * ( 0.01*inv_dx*(1+x).* solveFEM_assign3_phi_x_L( xi ) * solveFEM_assign3_phi_x_L( xi ) ...
                + 1.01 * solveFEM_assign3_phi_x_L( xi ) * solveFEM_assign3_phi_L( xi ) );
        if ( iR > 0 )
            % left to right
            T( iL, iR ) = T( iL, iR ) + w * ( 0.01*inv_dx*(1+x).* solveFEM_assign3_phi_x_R( xi ) * solveFEM_assign3_phi_x_L( xi ) ...
                    + 1.01 * solveFEM_assign3_phi_x_R( xi ) * solveFEM_assign3_phi_L( xi ) );
        end;
    end;
    
    if ( iR > 0 )
        % right to right
        T( iR, iR ) = T( iR, iR ) + w * ( 0.01*inv_dx*(1+x).* solveFEM_assign3_phi_x_R( xi ) * solveFEM_assign3_phi_x_R( xi ) ...
                + 1.01 * solveFEM_assign3_phi_x_R( xi ) * solveFEM_assign3_phi_R( xi ) );
        if ( iL > 0 )
            T( iR, iL ) = T( iR, iL ) + w * ( 0.01*inv_dx*(1+x).* solveFEM_assign3_phi_x_L( xi ) * solveFEM_assign3_phi_x_R( xi ) ...
                    + 1.01 * solveFEM_assign3_phi_x_L( xi ) * solveFEM_assign3_phi_R( xi ) );
        end;
    end;
    
end;

% compute the tight hadn side, we only have to add the boundary conditions
% for the right boundary
F = zeros( N-1, 1 );
[ w, x, xi ] = solveFEM_assign3_gauss3point( mesh_points(N,1), mesh_points(N,2) );
inv_dx = 2 / ( mesh_points(N,2) - mesh_points(N,1) );
F(N-1) = F(N-1) - w * ( 0.01*inv_dx*(1+x).* solveFEM_assign3_phi_x_R( xi ) * solveFEM_assign3_phi_x_L( xi ) ...
            + 1.01 * solveFEM_assign3_phi_x_R( xi ) * solveFEM_assign3_phi_L( xi ) );

% solve the system
C = T \ F;

% add back the boundary conditions for the solution
C = [ 0; C; 1; ];
%x = 0:1/N:1;

N_coarse = floor( N / 2 );
N_dense  = N - N_coarse;
x = [ 0:(break_point/N_coarse):break_point, break_point+((1-break_point)/N_dense):((1-break_point)/N_dense):1 ];

plot( x, C' );


function [ mesh_points, mesh_index ] = solveFEM_assign3_formUniformMesh( N )
% this creates a uniform mesh over the domain (0,1), the number of elements is N
% mesh_points is Nx2 numbers giving (a_i,b_i)
% mesh_index is Nx2 numbering specifying the indexes for the test funciton
% at each point of the mesh
x = 0:(1/N):1;

mesh_points = size( N, 2 );
mesh_index = size( N, 2 );

for k = 1:N
    mesh_points( k, 1 ) = x(k);
    mesh_points( k, 2 ) = x(k+1);
    
    mesh_index( k, 1 ) = k-1;
    mesh_index( k, 2 ) = k;
end;

mesh_index( 1, 1 ) = -1; % the first and last point will not
mesh_index( N, 2 ) = -1;


function [ mesh_points, mesh_index ] = solveFEM_assign3_formNonUniformMesh( N, break_point )
% this creates a uniform mesh over the domain (0,1), the number of elements is N
% mesh_points is Nx2 numbers giving (a_i,b_i)
% mesh_index is Nx2 numbering specifying the indexes for the test funciton
% at each point of the mesh
N_coarse = floor( N / 2 );
N_dense  = N - N_coarse;

x = [ 0:(break_point/N_coarse):break_point, break_point+((1-break_point)/N_dense):((1-break_point)/N_dense):1 ];

mesh_points = size( N, 2 );
mesh_index = size( N, 2 );

for k = 1:N
    mesh_points( k, 1 ) = x(k);
    mesh_points( k, 2 ) = x(k+1);
    
    mesh_index( k, 1 ) = k-1;
    mesh_index( k, 2 ) = k;
end;

mesh_index( 1, 1 ) = -1; % the first and last point will not
mesh_index( N, 2 ) = -1;

function [ w, x, xi ] = solveFEM_assign3_gauss3point( a, b )
w  = [ 5/9, 8/9, 5/9 ];
xi = [ -sqrt(3/5); 0; sqrt(3/5); ];
x = 0.5*(b+a) + 0.5*(b-a).*xi;

function [ y ] = solveFEM_assign3_phi_L( xi )
y = 0.5 - 0.5*xi;
function [ y ] = solveFEM_assign3_phi_R( xi )
y = 0.5 + 0.5*xi;
function [ y ] = solveFEM_assign3_phi_x_L( xi )
y = -0.5;
function [ y ] = solveFEM_assign3_phi_x_R( xi )
y = 0.5;
