Minimal surface finite differences problem - Matlab assemble

I face to the following problem:

$$(1+u_x^2)u_{yy} - 2u_xu_yu_{xy} + (1+u_y^2)u_{xx}=0.$$ Problem needs to be discretized and assembled. Does anybody know how to proceed in Matlab?

• Have you checked Concus, Paul. "Numerical solution of the minimal surface equation." Mathematics of Computation (1967): 340-350.? JSTOR link – nicoguaro Aug 19 '15 at 20:46
• What makes the problem more difficult in <put your language here>? – nicoguaro Aug 20 '15 at 19:24
• @nicoguaro in <Matlab> ... nothing more difficult - I guess it is same difficulty, I am just confused... after assemble, do I end up with one matrix and following system Ax = b? (I know - stupid question - I am really newbie in PDE and discretization in general)... also, afterwords - I need to obtain hessian and gradient - do u have any suggestions? – SmallElephant Aug 20 '15 at 19:41
• Hessian of gradient of what function? Have you solved any PDE using finite differences before? – nicoguaro Aug 20 '15 at 19:42
• seems like it is, you can write it in form as : $$A = \int \int (1 + u_x^2 + u_y^2)^{1/2} dx dy$$ and according to : en.wikipedia.org/wiki/Minimal_surface#cite_note-Meusnier1785-3 it should be – SmallElephant Aug 20 '15 at 20:17

Here is a fast implementation using sparse matrices and sparse Jacobian estimation:

% 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);
% use Kronecker product to build matrix of d/dx and d/dy
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",...