function x = solve_ldl(A,b)

% Compute the LDL^T factorization of the matrix A using the program
% myldl.m.
LD = myldl(A);

n = length(A);

% initialize the solution x as zero vector of the correct size
x = zeros(size(b));

% Solve the equation Ly = b, call the solution again b.
for k=1:n-1
    b(k+1:end) = b(k+1:end)-LD(k+1:end,k)*b(k);
end

% Divide the result componentwise by the diagonal entries of LD.
b = b./diag(LD);

% Solve the equation L^T x = b.
for i=n:-1:1
    x(i) = b(i)-sum(LD(i+1:end,i).*x(i+1:end));
end