% Finite element script to solve the heat equation
% d/dt u - Lap(u) = f
% on the unit circle.

clear all; clc; close all;
profile on

addpath('distmesh') % adds the toolbox for space discretization

%% Initializations
hs = [0.2 0.1 0.07 0.05 0.03 0.02 0.01];
hs_computed = zeros(length(hs),1);
taus = [0.5 0.25,0.2,0.125,0.1,0.05, 0.02];
T_end = 4;
errors_L2_BDF = zeros(length(hs),length(taus));
errors_H1_BDF = zeros(length(hs),length(taus));
errors_L2_RadauII = zeros(length(hs),length(taus));
errors_H1_RadauII = zeros(length(hs),length(taus));
fd=@(p) sqrt(sum(p.^2,2))-1; % distance function for the grid generation

%% Iterate over all grid sizes h
for h_count = 1:length(hs)
    h = hs(h_count);
    %% Create grid
    tic
    [nodes,elements]=distmesh2d(fd,@huniform,hs(h_count),[-1,-1;1,1],[]);
    time = toc;
    fprintf('Mesh generation for h = %f completed after %f seconds.\n',hs(h_count),time);
    fprintf('The mesh has %d nodes and %d elements.\n',size(nodes,1),size(elements,1));
    
    %% Compute the 'real' grid size
    % ... compute the real grid size of the grid generated by distmesh and
    % store it in the vector hs_computed
    
    %% Determine boundary points and interior points
    % determine indices of the boundary points and interior points and
    % store them in vectors ...
    
    %% Assembly mass and stiffness matrix
    % use your own matrix assembly code from numerics 3. If you don't have
    % it, you can use the given function.
   
    %% Iterate over all tau
    for tau_count = 1:length(taus) % iterate over all tau
        %% Compute 2 initial values using the analytical solution
        % ...
        
        %% BDF2 method
        % ... BDF 2 method goes here.
        
        %% RadauII method
        % ... RadauII method goes here.
        
        %% Compute error
        % .. compute errors of the numerical solution and store them in the
        % error matrices
    end
end

%% Generate convergence plots
% generate 4 convergence plots per method in one window using subplot. 
% Remind to use the real grid size instead of the distmesh input parameters.

profile off
profile viewer