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}];