What' s the problem in the following code? Could you help me correct the code?

restart:
with(LinearAlgebra):
with(plots):
k:=2:M:=2:
Tm:=proc(t,m)
option remember;
if m=0 then return 1 else
if m=1 then return t else
2*t*Tm(t,m-1)-Tm(t,m-2) fi; fi;
end proc:
alpha:=(m)->piecewise(m=0,1/sqrt(Pi),sqrt(2)/sqrt(Pi)):
Phi:=(n,m,x)->piecewise((n-1)/(2^(k-1)) <= x and x <= n/(2^(k-1)),
alpha(m)*(2^(k/2))*Tm2(2^k*x-2*n+1,m), 0):
Phi2:=(t)->Array([seq(seq(Phi(i,j,t),j=0..M-1),i=1...2^(k-1))] ):
Phi2 := t -> Array([seq(seq(Phi(i,j,t),j = 0 .. M-1),i = 1 .. 2^(k-1))]):
for i from 1 to ((2^(k-1))*M) do
r[i]:=evalf(Phi2((2*i-1)/((2^k)*M)));
end do:
Phi_mxm:=Matrix([seq([r[i]],i=1...(2^(k-1))*M)])^+:
Phi_mxm_Inv:=MatrixInverse(Phi_mxm): m:=M*(2^(k-1)):
PB_nn:=Matrix(m,m):
i:=(i,n)->((i+1)^(n+1)-2*i^(n+1)+(i-1)^(n+1)):
nn:=1:#xi indisi
for i from 1 to m do p:=0:
for j from 1 to m do
if i=j then PB_nn[i,j]:=1:
fi:
if i<j then p:=p+1:
PB_nn[i,j]:= xi(p,nn):
fi:
end do:
p:=1:
end do:
for i from 1 to m do p:=0:
for j from 1 to m do
PB_nn:=(1/(m^nn))*(1/((nn+1)!))* PB_nn:
P_nn1:=Phi_mxm.PB_nn.MatrixInverse(Phi_mxm):
if i<j then p:=p+1:
PB_nn[i,j]:= xi(p,nn):
fi:
end do:
p:=1:
end do:
PB_nn:=(1/(m^nn))*(1/((nn+1)!))* PB_nn:
P_nn2:=Phi_mxm.PB_nn.MatrixInverse(Phi_mxm):
nn:=3:#xi indisi
for i from 1 to m do p:=0:
for j from 1 to m do
if i=j then
PB_nn[i,j]:=1: fi: if i<j then p:=p+1:
B_nn[i,j]:= xi(p,nn):
fi:
end do:
p:=1:
end do:
PB_nn:=(1/(m^nn))*(1/((nn+1)!))* PB_nn:
for j from 1 to (2^(k-1))*M do
x[j]:=evalf((2*j-1)/((2^k)*M)):
end:
C:=Matrix( (2^(k-1))*M,(2^(k-1))*M,symbol=c):
for i from 1 to (2^(k-1))*M do
for j from 1 to (2^(k-1))*M do
PDE(i,j):=simplify(evalf(Phi2(x[i]).(P_nn3).C.(Phi2(t[j])^+)+x[i]/(6*((1-t[j])^2))-x[i]*Phi2(1).(P_nn3).C.(Phi2(t[j])^+)-(1/(6*((1-t[j])^2)))-6*(Phi2(x[i]).(P_nn3).C.(Phi2(t[j])^+)+x[i]/(6*((1-t[j])))-x[i]*Phi2(1).(P_nn3).C.(Phi2(t[j])^+)-(1/(6*((1-t[j]))))).(Phi2(x[i]).(P_nn2).C.(Phi2(t[j])^+)+(1/(6*(1-t[j])))-Phi2(1).(P_nn3).C.(Phi2(t[j])^+))+Phi2(x[i]).C.(P_nn1).(Phi2(t[j])^+)))=0 end do;
end do;
sys:=[seq(seq(PDE(i,j),i=1..(2^(k-1))*M ),j=1..(2^(k-1))*M )]:
var:=[seq(seq(c[i,j],j=1..(2^(k-1))*M ),i=1..(2^(k-1))*M )]:
A,b:=GenerateMatrix(sys,var);
G:=Matrix((2^(k-1))*M*(2^(k-1))*M,1,symbol=g);
A.G=b;
coz:=fsolve({seq(seq(PDE(i,j),i=1..(2^(k-1))*M ),j=1..(2^(k-1))*M )}):
C:=subs(op(coz),C);
u:=(x,t)->(1/6)*((x-1)/(1-t)):
U_ERR:=(x,t)->abs(U(x,t)-u(x,t)):
xbas:=1:
tbas:=1:
with(plots):
g1:=plot3d(U(x,t),x=0..xbas,t=0..tbas,color=red):
g2:=plot3d(u(x,t),x=0..xbas,t=0..tbas):
display(g1,g2);