Question: plotting error

Any one help me  to remove the error.

I want to plot the curve for different values of alpha.

here is my codes.

thanks in advance 

 

restart:
with(linalg):with(plots): 
ge[1]:=diff(u[1](x,t),t)=alpha*diff((u[2](x,t)-1)*diff(u[1](x,t),x),x)+(16*x*t-2*t-16*(u[2](x,t)-1))*(u[1](x,t)-1)+10*x*exp(-4*x):
ge[2]:=diff(u[2](x,t),t)=diff(u[2](x,t),x$2)+alpha*diff(u[1](x,t),x)+4*(u[1](x,t)-1)+x^2-2*t-10*t*exp(-4*x): 
bc1[1]:=u[1](x,t)-1: 
bc1[2]:=u[2](x,t)-1: 
bc2[1]:=3*u[1](x,t)+diff(u[1](x,t),x)-3:
bc2[2]:=5*diff(u[2](x,t),x)-evalf(exp(4))*(u[1](x,t)-1):
IC[1]:=u[1](x,0)=1: 
IC[2]:=u[2](x,0)=1: 
NN:=2: 
N:=2:
L:=1:
for i to NN do  
dydxf[i]:=1/2*(-u[2,i](t)-3*u[0,i](t)+4*u[1,i](t))/h: 
dydxb[i]:=1/2*(u[N-1,i](t)+3*u[N+1,i](t)-4*u[N,i](t))/h:
dydx[i]:=1/2/h*(u[m+1,i](t)-u[m-1,i](t)); 
d2ydx2[i]:=1/h^2*(u[m-1,i](t)-2*u[m,i](t)+u[m+1,i](t)):od:
 for i to NN do bc1[i]:=subs(diff(u[1](x,t),x)=dydxf[1],
diff(u[2](x,t),x)=dydxf[2],u[1](x,t) 
=u[0,1](t),u[2](x,t)=u[0,2](t),x=0,bc1[i]):od: 
for i to NN do bc2[i]:=subs(diff(u[1](x,t),x)=dydxb[1],
diff(u[2](x,t),x)=dydxb[2],u[1](x,t) 
=u[N+1,1](t),u[2](x,t)=u[N+1,2](t),x=L,bc2[i]):od:
for i to NN do eq[0,i]:=bc1[i];eq[N+1,i]:=bc2[i]:od: 
for i from 1 to N do eq[i,1]:=diff(u[i,1](t),t)= subs(diff(u[1](x,t),x$2) =
subs(m=i,d2ydx2[1]), 
diff(u[2](x,t),x$2) = subs(m=i,d2ydx2[2]),diff(u[1](x,t),x) =
subs(m=i,dydx[1]),diff(u[2](x,t),x) = subs(m=i,dydx[2]),u[1](x,t)=u[i,1](t), 
u[2](x,t)=u[i,2](t),x=i*h,rhs(ge[1])):od:

for i from 1 to N do eq[i,2]:=diff(u[i,2](t),t)= subs(diff(u[1](x,t),x$2) =
subs(m=i,d2ydx2[1]), 
diff(u[2](x,t),x$2) = subs(m=i,d2ydx2[2]),diff(u[1](x,t),x) =
subs(m=i,dydx[1]),diff(u[2](x,t),x) = subs(m=i,dydx[2]),u[1](x,t)=u[i,1](t),
u[2](x,t)=u[i,2](t),x=i*h,rhs(ge[2])):od: 

for i to NN do u[0,i](t):=(solve(eq[0,i],u[0,i](t))):od:

 for i to NN do u[N+1,i](t):=(solve(eq[N+1,i],u[N+1,i](t))):od:

 h:=L/(N+1): 

for i from 1 to N do eq[i,1]:=eval(eq[i,1]):od: 
 for i from 1 to N do eq[i,2]:=eval(eq[i,2]):od:

eqs:=seq(seq((eq[i,j]),i=1..N),j=1..NN): 
Y:=seq(seq(u[i,j](t),i=1..N),j=1..NN): 

 ICs:=seq(u[i,1](0)=rhs(IC[1]),i=1..N),seq(u[i,2](0)=rhs(IC[2]),i=1..N): 

sol:=dsolve({eqs,ICs},{Y},type=numeric,stiff=true,maxfun=1000000,abserr=1e-6,relerr=1e-5,output=listprocedure):

Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)
for j to NN do for i to N do U[i,j]:=subs(sol,u[i,j](t)):od:od: 

for i to NN do U[0,i]:=subs(u[1,1](t)=U[1,1],u[1,2](t)=U[1,2],
u[2,1](t)=U[2,1],u[2,2](t)=U[2,2],u[0,i](t)):od:
 for i to NN do U[N+1,i]:=eval(subs(u[N,1](t)=U[N,1],u[N,2](t)=U[N,2],
u[N-1,1](t)=U[N-1,1],u[N-1,2](t)=U[N-1,2],u[N+1,i](t))):od:
tf:=1.: 
M:=30: 
T1:=[seq(tf*i/M,i=0..M)]: 
PP:=matrix(N+2,M+1): 
for i from 1 to N+2 do PP[i,1]:=evalf(subs(x=(i-1)*h,rhs(IC[1]))):od: 
for i from 1 to N+2 do for j from 2 to M+1 do
PP[i,j]:=evalf(subs(t=T1[j],U[i-1,1](t))):od:od:
 
G1:=[seq([ seq([(i-1)*h,T1[j],PP[i,j]], i=1..N+2)], j=1..M+1)]: 
t=0.02: 
pars:=[0.1,0.5,1,2,5];
clr:=[black,red,green,gold,blue];  
for m from 1 to 5 do 
G1[m]:=plot([seq(subs(alpha=pars[m],G1[i])],i=0..N+1)],thickness=3,color=clr[j]):od:  

display({seq(G1[i],i=1..5)},title="Figure ",axes=boxed,labels=[x,u]);
restart:
 

Please Wait...