  function ShowLDLIndef()
% Using ldl to solve a symmetric indefinite system...
clc
format short
n = 7;
A = randn(n,n); A = A + A';
b = randn(n,1);
[L,D,P] = ldl(A)

% P'AP = LDL', i.e., A = PLDL'P'
% To solve Ax = b  we want to solve P(L(D(L'(P'x)))) = b

x = P*(L'\(D\(L\(P'*b))));

% D is now block diagonal with 1x1 and 2x2 diagonal blocks