# Orthogonal Collocation Using Jacobi Polynomials

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 = 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;
Table[C[i] == 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'[t] == 0,C[t] == 1, Table[C[i] == 10, {i, 2, 9}]}],Table[C[i][t], {i, 1, 10}],{t, 0, 10}];

• I don't currently have Mathematica at hand, so can't run your code but I can see a couple of potential issues: you use the reserved variable C (better not use variables starting with capitals at all). Furthermore, it looks like the three Table lines at the bottom don't do anything at all, at least the result isn't assigned to a variable. You could also try to change the 1.0 to 1 in the second line. – Sjoerd C. de Vries May 11 '15 at 11:25
• Please show what's returned in sol at the end. – Sjoerd C. de Vries May 11 '15 at 11:35
• I corrected the use of C, a reserved variable, and removed the redundant tables. I still get the same issue with the ndsolve. Mathematica says: Equation or list of equations expected instead of... – ABDUL-RASHID BAWAH May 12 '15 at 16:41

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[C[i]][t] == piece1[[i]] + piece2[[i-1]], {i, 2, 9}],
Derivative[C][t] == 0,
Derivative[C][t] == 0,
(*I fixed these two constraints.
Not sure what they are meant to be, change if necessary*)
C == 1,
C == 1,
Table[C[i] == 10, {i, 2, 9}]}];

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