Orthogonal Collocation Using Jacobi Polynomials

Question

I'm trying to solve a PDE(diffusion-reaction in a spherical catalyst pellet) using Jacobi Orthogonal Collocation method. But at the stage of solving the system of ODEs(using NDSolve) resulting from the collocation, Mathematica does not return a result. This is the equation I am trying to solve: $$\frac{\partial C(t,r)}{\partial t} = \frac{\partial ^2C(t,r)}{\partial r^{2}}+\frac{1}{r}\frac{\partial C(t,r)}{\partial r}+\frac{C(t, r)}{(C(t,r)+1)^2}$$ You can find the code I have so far below:

Np = 10;
sol = Solve[JacobiP[Np - 2, 0, 0, 2 x - 1.] == 0, x]
tbl = Select[Table[x /. sol[[i]], {i, 1, Length[sol]}], # > 0 &];
x[Np] = 1.0;
x[1] = 10^-20; For[i = 2, i < Np, {x[i] = tbl[[i - 1]], i++}];
Q = Chop[Table[Table[x[j]^(i - 1), {i, 1, Np}], {j, 1, Np}]];
IQ = Inverse[Q] // Quiet;
Ci = Chop[Table[Table[(i - 1) x[j]^(i - 2), {i, 1, Np}], {j, 1, Np}]];
A = Ci.IQ;
Di = Chop[
Table[Table[(i - 1) (i - 2) x[j]^(i - 3), {i, 1, Np}], {j, 1, 
 Np}]];
    B = Di.IQ;
B // MatrixForm;
Z = IdentityMatrix[8];
Table[C[i][0] == 10, {i, 2, 9}]
Table[C'[i][t] == (B.Table[C[i][t], {i, 1, 10}])[[i]], {i, 2, 9}] +Z.Table[C[i][t]/(1 + C[i][t])^2, {i, 2, 9}];
 Flatten[Table[C'[i][t] == (B.Table[C[i][t], {i, 1, 10}])[[i]], {i, 2, 9}] +Z.Table[C[i][t]/(1 + C[i][t])^2, {i, 2, 9}]];
sol = NDSolve[Flatten[{Table[C'[i][t] == (B.Table[C[i][t], {i, 1, 10}])[[i]], {i, 2, 9}] +Z.Table[C[i][t]/(1 + C[i][t])^2, {i, 2, 9}], C'[0][t] == 0,C[1][t] == 1, Table[C[i][0] == 10, {i, 2, 9}]}],Table[C[i][t], {i, 1, 10}],{t, 0, 10}];

Answer 1

It's being caused by a misplaced ')'. Nested tables are frequently a bad idea. Here:

Table[C'[i][t] == (B.Table[C[i][t], {i, 1, 10}])[[i]], {i, 2, 9}] + ...

You can't do that with the plus on the end. In the equation you're feeding into NDSolve, essentially you have (stuff == moreStuff) + extra which makes no sense. Once I correct for that the error goes away. Manually inspect your equations. I fixed up some of your constraints as the system was over/underdetermined so you may want to check this more closely:

vars = Table[C[i][t], {i, 0, 9}];
piece1 = B.vars[[1 ;; 10]];
piece2 = Z.Table[C[i][t]/(1 + C[i][t])^2, {i, 2, 9}];
equations = 
  Flatten[{Table[
   Derivative[1][C[i]][t] == piece1[[i]] + piece2[[i-1]], {i, 2, 9}],
   Derivative[1][C[1]][t] == 0,
   Derivative[1][C[0]][t] == 0,
   (*I fixed these two constraints. 
    Not sure what they are meant to be, change if necessary*)
   C[0][0] == 1,
   C[1][0] == 1,
   Table[C[i][0] == 10, {i, 2, 9}]}];

sol = NDSolve[equations, vars, {t, 0, 10}];