function [a,o,n,b,x] = substitute(a,o,n,b);
%first foward substitute
%solve [L]{d}={b} 
%note that d is compactly stored in {b} as a result
%we start from i=2 and note that d(1) is b(1)


for i = 2 : n
    sum=b(o(i));
    for j=1: i-1
        sum=sum - a(o(i),j)*b(o(j));
    end
    b(o(i))=sum;
end
b
x(n)=b(o(n))/a(o(n),n);

%now back wards substitute
%solve [U]{x}={d}

for i = n-1 :-1 :1
    sum=0;
    for j = i+1:n
        sum=sum+a(o(i),j)*x(j);
    end
    x(i)=(b(o(i))-sum)/a(o(i),i);
end

b