Dear Maple researchers

I have a problem in solving a system of odes that resulted from discretizing, in space variable, method of lines (MOL).

The basic idea of this code is constructed from the following paper:

http://www.sciencedirect.com/science/article/pii/S0096300313008060

If kindly is possible, please tell me whas the solution of this problem.

With kin dregards,

Emran Tohidi.

My codes is here:

> restart;

> with(orthopoly);

print(`output redirected...`); # input placeholder

> N := 4; Digits := 20;

print(`output redirected...`); # input placeholder

> A := -1; B := 1; rho := 3/4;

> g1 := proc (t) options operator, arrow; 1/2+(1/2)*tanh((1/2)*(A-(2*rho-1)*t/sqrt(2))/sqrt(2)) end proc; g2 := proc (t) options operator, arrow; 1/2+(1/2)*tanh((1/2)*(B-(2*rho-1)*t/sqrt(2))/sqrt(2)) end proc;

print(`output redirected...`); # input placeholder

> f := proc (x) options operator, arrow; 1/2+(1/2)*tanh((1/2)*x/sqrt(2)) end proc;

print(`output redirected...`); # input placeholder

> uexact := proc (x, t) options operator, arrow; 1/2+(1/2)*tanh((1/2)*(x-(2*rho-1)*t/sqrt(2))/sqrt(2)) end proc;

print(`output redirected...`); # input placeholder

> basiceq := simplify(diff(uexact(x, t), `$`(t, 1))-(diff(uexact(x, t), `$`(x, 2)))+uexact(x, t)*(1-uexact(x, t))*(rho-uexact(x, t)));

print(`output redirected...`); # input placeholder

0

> alpha := 0; beta := 0; pol := P(N-1, alpha+1, beta+1, x); pol := unapply(pol, x); dpol := simplify(diff(pol(x), x)); dpol := unapply(dpol, x);

print(`output redirected...`); # input placeholder

> nodes := fsolve(P(N-1, alpha+1, beta+1, x));

%;

> xx[0] := -1;

> for i to N-1 do xx[i] := nodes[i] end do;

print(`output redirected...`); # input placeholder

> xx[N] := 1;

> for k from 0 to N do h[k] := 2^(alpha+beta+1)*GAMMA(k+alpha+1)*GAMMA(k+beta+1)/((2*k+alpha+beta+1)*GAMMA(k+1)*GAMMA(k+alpha+beta+1)) end do;

print(`output redirected...`); # input placeholder

> w[0] := 2^(alpha+beta+1)*(beta+1)*GAMMA(beta+1)^2*GAMMA(N)*GAMMA(N+alpha+1)/(GAMMA(N+beta+1)*GAMMA(N+alpha+beta+2));

print(`output redirected...`); # input placeholder

> for jj to N-1 do w[jj] := 2^(alpha+beta+3)*GAMMA(N+alpha+1)*GAMMA(N+beta+1)/((1-xx[jj]^2)^2*dpol(xx[jj])^2*factorial(N-1)*GAMMA(N+alpha+beta+2)) end do;

print(`output redirected...`); # input placeholder

> w[N] := 2^(alpha+beta+1)*(alpha+1)*GAMMA(alpha+1)^2*GAMMA(N)*GAMMA(N+beta+1)/(GAMMA(N+alpha+1)*GAMMA(N+alpha+beta+2));

print(`output redirected...`); # input placeholder

> for j from 0 to N do dpoly1[j] := simplify(diff(P(j, alpha, beta, x), `$`(x, 1))); dpoly1[j] := unapply(dpoly1[j], x); dpoly2[j] := simplify(diff(P(j, alpha, beta, x), `$`(x, 2))); dpoly2[j] := unapply(dpoly2[j], x) end do;

print(`output redirected...`); # input placeholder

print(??); # input placeholder

> for n to N-1 do for i from 0 to N do BB[n, i] := sum(P(jjj, alpha, beta, xx[jjj])*dpoly2[jjj](xx[n])*w[i]/h[jjj], jjj = 0 .. N) end do end do;

> for n to N-1 do d[n] := BB[n, 0]*g1(t)+BB[n, N]*g2(t); d[n] := unapply(d[n], t) end do;

print(`output redirected...`); # input placeholder

> for nn to N-1 do F[nn] := simplify(sum(BB[nn, ii]*u[ii](t), ii = 1 .. N-1)+u[nn](t)*(1-u[nn](t))*(rho-u[nn](t))+d[nn](t)); F[nn] := unapply(F[nn], t) end do;

print(`output redirected...`); # input placeholder

> sys1 := [seq(d*u[q](t)/dt = F[q](t), q = 1 .. N-1)];

print(`output redirected...`); # input placeholder

[d u[1](t)

[--------- = 40.708333333333333334 u[1](t) + 52.190476190476190476 u[2](t)

[ dt

2 3

+ 39.958333333333333334 u[3](t) - 1.7500000000000000000 u[1](t) + u[1](t)

+ 7.3392857142857142858

- 3.6696428571428571429 tanh(0.35355339059327376220

+ 0.12500000000000000000 t) - 3.6696428571428571429 tanh(

d u[2](t)

-0.35355339059327376220 + 0.12500000000000000000 t), --------- =

dt

-20.416666666666666667 u[1](t) - 25.916666666666666667 u[2](t)

2 3

- 20.416666666666666667 u[3](t) - 1.7500000000000000000 u[2](t) + u[2](t)

- 3.7500000000000000000

+ 1.8750000000000000000 tanh(0.35355339059327376220

+ 0.12500000000000000000 t) + 1.8750000000000000000 tanh(

d u[3](t)

-0.35355339059327376220 + 0.12500000000000000000 t), --------- = 29.458333333\

dt

333333333 u[1](t) + 38.476190476190476190 u[2](t)

2 3

+ 30.208333333333333333 u[3](t) - 1.7500000000000000000 u[3](t) + u[3](t)

+ 5.4107142857142857144

- 2.7053571428571428572 tanh(0.35355339059327376220

+ 0.12500000000000000000 t) - 2.7053571428571428572 tanh(

]

-0.35355339059327376220 + 0.12500000000000000000 t)]

]

> ics := seq(u[qq](0) = evalf(f(xx[qq])), qq = 1 .. N-1);

print(`output redirected...`); # input placeholder

u[1](0) = 0.38629570659055483825, u[2](0) = 0.50000000000000000000,

u[3](0) = 0.61370429340944516175

> dsolve([sys1, ics], numeic);

%;

Error, (in dsolve) invalid input: `PDEtools/sdsolve` expects its 1st argument, SYS, to be of type {set({`<>`, `=`, algebraic}), list({`<>`, `=`, algebraic})}, but received [[d*u[1](t)/dt = (20354166666666666667/500000000000000000)*u[1](t)+(13047619047619047619/250000000000000000)*u[2](t)+(19979166666666666667/500000000000000000)*u[3](t)-(7/4)*u[1](t)^2+u[1](t)^3+36696428571428571429/5000000000000000000-(36696428571428571429/10000000000000000000)*tanh(1767766952966368811/5000000000000000000+(1/8)*t)-(36696428571428571429/10000000000000000000)*tanh(-1767766952966368811/5000000000000000000+(1/8)*t), d*u[2](t)/dt = -(20416666666666666667/1000000...