这是使用稀疏矩阵和稀疏雅可比估计的快速实现：

% define square domain [-1,1] x [-1,1]
n = 51;
x=linspace(-1,1,n);
y=x;
[X,Y]=meshgrid(x,x);

% build finite differences operators
dx=x(2)-x(1);
e=ones(n,1);
d0x=ones(n,1);
grad = spdiags([-e 0*e e],-1:1,n,n)/2/dx;
% use Kronecker product to build matrix of d/dx and d/dy
gradx = kron(grad,speye(n,n));
grady = kron(speye(n,n),grad);
lap = spdiags([e -2*e e],-1:1,n,n)/dx^2;
% use Kronecker product to build matrix of d/dx^2 and d/dy^2
lapx = kron(lap,speye(n,n));
lapy = kron(speye(n,n),lap);

% Dirichlet boundary condition
bdy = find(X(:)==x(1) | X(:)==x(end) | Y(:)==y(1) | Y(:)==y(end));
cnd = 0.5*sign(sin(4*(X(bdy)+Y(bdy))));

% initial value of iterate
f0 = zeros(n^2,1)

opt = optimoptions("fsolve","Algorithm","trust-region",...
    "JacobPattern",lapy+lapx+grady*gradx,"display","iter");

f = fsolve(@(f) fun(f,gradx,grady,lapx,lapy,bdy,cnd),f0,opt);

surf(x,x,reshape(f,n,n),"facecolor","interp")

function out=fun(f,gradx,grady,lapx,lapy,bdy,cnd)
    fx = gradx*f;
    fy = grady*f;
    % equation of mininal surface. Surface is the graph of z=f(x,y)
    % inside domain
    out = (1+fy.^2).*(lapx*f)+(1+fx.^2).*(lapy*f)-2*fx.*fy.*(gradx*(grady*f));
    % at the boundary
    out(bdy) = f(bdy)-cnd;
end