makeGram:=proc(T,lvl,funclist)
	
	local G, size, i, j;
	
	size:=nops(funclist);
	G:=matrix(size,size);

	for i from 1 to size do
		for j from 1 to size do
				
		G[i,j]:=discInnerProduct(T,lvl,funclist[i],funclist[j]);

		end do;
	end do;

	return(G);

end;

makeDiagInside:=proc(G)
	local D, i,e,eig,l;

	#D:=diag(evalf(eigenvals(G)));
	l:=[];
	eig:=evalf(eigenvects(G));
	for e in eig do
                for i from 1 to eig[2] do
              		l:=[op(l),eig[1]];          
                end do;
        end do;

	D:=diag(l);

	for i in indices(D) do
		D[op(i)]:=1/sqrt(D[op(i)]);
	end do;
	return(D);
end;

makeDiagOutside:=proc(G)
	local out, col, column, eig, e, i;

	column:=0;
	col:=[];
	
	eig:=evalf(eigenvects(G));

	for e in eig do
		for i from 1 to e[2] do	
			col:=[op(col),e[3][i]];
			column:=column + 1;
		end do;	
	end do;
	
	out:=matrix(column, column, col);
	out:=transpose(out);
	return(out);

end;

makeA:=proc(T,lvl,funclist)
	local U, D, G, A;

	G:=makeGram(T,lvl,funclist);
	U:=makeDiagOutside(G);
	D:=makeDiagInside(G);
	A:=multiply(U,D,inverse(U));

	return A;

end;

makeOrtho:=proc(T, lvl, funclist, startfunc)
	local A, i, j;

	A:=makeA(T, lvl, funclist);

	for i from 1 to nops(funclist) do
		clearFunc(T,lvl,startfunc + i -1); 
		for j from 1 to nops(funclist) do
			funcSumWithConstant(T,lvl,A[i,j],startfunc + i -1,funclist[j] ,startfunc + i -1);
		end do;
	end do;

	return(startfunc + nops(funclist));
end;