
function femplotexample

    %% dimensions of domain
    Lx = 10.0;
    Ly = 3.0;

    % discretisation - number of elements in x,y
    mx = 4;
    my = 5;

    nx = mx + 1;
    ny = my + 1;

    dx = Lx/mx;
    dy = Ly/my;


    %% assign mesh coordinates
    gcoords = zeros(2,nx*ny);

    for i = 1: nx
        for j = 1 : ny
            xn = (i-1)*dx;
            yn = (j-1)*dy;
            nid = (i-1)+(j-1)*nx + 1;
            gcoords(1,nid) = xn;
            gcoords(2,nid) = yn;
        end
    end

    %% pertub on node at the top right side

    gcoords(2,nx*ny-2) = Ly + 0.1;


    %% solution vector
    phi = zeros(1,nx*ny);

    %% initial condition
    for i = 1 : nx
        for j = 1 : ny
            nid = (i-1)+(j-1)*nx+1;

            xn = gcoords(1,nid);
            yn = gcoords(2,nid);
            phi(nid) = exp(-( (xn-2.5)^2 + (yn-3.0)^2 ) );
        end
    end

    plot_fem_solution(nx,ny,gcoords,phi);

end


%% --------------------------------------------------------------
%% basic plotting function
    function plot_fem_solution(nx,ny,coordinates,phi)			
    %% transform coords into mesh grid format
    Xc = zeros(nx,ny);
    Yc = zeros(nx,ny);
    %% transform data into mesh grid format							 
    phic = zeros(nx,ny);
    for i = 1: nx
        for j = 1 : ny
            nid = (i-1)+(j-1)*nx + 1;
            
            Xc(i,j) = coordinates(1,nid);
            Yc(i,j) = coordinates(2,nid);

            phic(i,j) = phi(nid);
        end
    end

    figure
    surf(Xc,Yc,phic);
end

