function minsurf
    xmin= 0; xmax= 1;
    ymin= 0; ymax= 1;
    n= 20; m= 20;
    hx=(xmax-xmin)/n; hy= (ymax-ymin)/m;
    xvalues= xmin:hx:xmax;
    yvalues= ymin:hy:ymax;
    numpoints= (m-1)*(n-1);
    p= zeros(1,numpoints); % initial values
    
    p= fminunc(@fdint,p);
    
    plots(p)
    function s=fdint(p)
        [x,y,z]= points_to_matrix(p);
        s= 0;
        for i=1:n
            for j=1:m
                dzdx= (z(i+1,j)-z(i,j))/hx;
                dzdy= (z(i,j+1)-z(i,j))/hy;
                s=s+hx*hy*sqrt(1+dzdx^2+dzdy^2);
            end
        end
    end
    function [x,y,z]= points_to_matrix(p)
        [x,y]= meshgrid(xvalues,yvalues);
        z= zeros(size(x));
        z(1,:)= zxmin(yvalues);
        z(n+1,:)= zxmax(yvalues);
        z(:,1)= zymin(xvalues);
        z(:,m+1)= zymax(xvalues);
        z(2:n,2:m)= reshape(p,[n-1 m-1]);
    end
    function plots(p)
        [x,y,z]= points_to_matrix(p);
        surf(x,y,z);
    end
end

function z= zxmax(y)
    z= 1+0*y;
end
function z= zxmin(y)
    z= 0*y;
end
function z= zymax(x)
    z= x.^2;
end
function z= zymin(x)
    z= x.^6;
end