% Eliminating process for column p
for k=p+1:N
m=Aug(k,p)/Aug(p,p);
Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);
end
end

%Back substitution on [U|Y]
X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));

The above program is to construct the solution to AX=B by reducing the augmented matrix [A B] to upper triangularization method and then performing back substitution.

I am trying to modify this program so that it will solve M linear systems with the same coefficient matrix A but different column matrices B . The M linear systems look like AX1 = B1 AX2 = B2 ...... AXm=Bm
(1 , 2 and m are in subscript)

% Eliminating process for column p
for k=p+1:N
m=Aug(k,p)/Aug(p,p);
Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);
end
end

%Back substitution on [U|Y]
X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));

The above program is to construct the solution to AX=B by reducing the augmented matrix [A B] to upper triangularization method and then performing back substitution.

I am trying to modify this program so that it will solve M linear systems with the same coefficient matrix A but different column matrices B . The M linear systems look like AX1 = B1 AX2 = B2 ...... AXm=Bm
(1 , 2 and m are in subscript)

Please help!!!

My guess is it should work as it is, just make B the matrix of the column vectors, all right you will have to modify the back substitution to do back sub for all the columns of augmentation (this can be done by looping over the extra columns, though the back substitution function might just work if you give it all the augmentation cols at once).