function []=muster();

format long;

clear all;
close all;

kmax=7;

% analytic solution u

u_ana=@(x,y)((x.^4)' * y.^5-17*sin(x'*y));

% define input functions f and g

fun=@(x,y)(-12* (x.^2)' * y.^5-20*(x.^4)' * y.^3-17*sin(x'*y) .* (repmat((x.^2)',1,length(y))+repmat(y.^2,length(x),1)));

gx0=@(y)(0);

gx1=@(y)(y.^5-17*sin(y));

gy0=@(x)(0);

gy1=@(x)(x.^4-17*sin(x));

% allocate space
N=zeros(1,kmax);
errsup=zeros(1,kmax);
err2=zeros(1,kmax);

% calculate discrete solutions and errors for different number of grid nodes

for k=1:kmax;
    N(k)=2^(k+1);
    % get analytic solution
    x=linspace(0,1,N(k)+1);
    y=linspace(0,1,N(k)+1);
    u_ana_mat=u_ana(x(2:N(k)),y(2:N(k)));
    u_ana_vec=u_ana_mat(:);
    
    % get discrete solution
    u=fdPoisson2d(N(k),fun,gx0,gy0,gx1,gy1);
    
    % calculate errors
    err_vec=u-u_ana_vec;
    
    errsup(k)=max(abs(err_vec));
    err2(k)= norm(err_vec);
end

% calculate and display convergence rates
rate_sup = (log(errsup(end)) - log(errsup(end-1)) ) / (log(1/N(end))-log(1/N(end-1)));
rate_2 = (log(err2(end)) - log(err2(end-1)) ) / (log(1/N(end))-log(1/N(end-1)));

s = sprintf('The order of convergence in the sup-norm based on finest meshes is %g',rate_sup);
disp(s);
t = sprintf('The order of convergence in the Euclidean norm based on finest meshes is %g',rate_2);
disp(t);

% plot errors

figure(2);
loglog(1./N,errsup, 1./N, err2,1./N, 1./N.^2, 1./N, 1./N);
legend('sup', 'l2', 'h^2', 'h');
ylabel('error');
xlabel('h');
title('logarithmic plot of the error dependent on the mesh width');

figure(3);
title('analytic solution on finest mesh')
surf(x,y,u_ana(x,y));