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: