How to solve numerically system of linear coupled PDE using Finite Difference?

Question

How to solve this coupled linear equation with the 2D plot of the variation of parameter with constant with I.C: u(y,0)=c(y,0)=theta(y,0)=0, B.C: u(0,t)=1,c(0,t)=theta(0,t)=t.

Answer 1

Since this forum is for using Mathematica to solve problems, here is how you would solve these using NDSolve.

If you want code using finite difference, I suggest the forum https://scicomp.stackexchange.com/ instead.

I fixed few things. First your BC and IC is not consistent, You say u(y,0)=0 but then say u(0,t)=1. So I changed u(0,t)=0.

You also do not have enough spatial B.C. You have derivatives w.r.t. y twice, but only provided one BC. So I added more BC's to make NDSolve happy. You can change these as needed.

You also did not provide numerical values for Gr,Gm,RmPr,Sc do I added dummy values for these. You can change them to the actual values.

ClearAll[u, y, t, c, θ, Gr, Gm, R, Pr, Sc]
pde1 = D[u[y, t], t] == Gr* θ[y, t] + Gm*c[y, t] + D[u[y, t], {y, 2}]
pde2 = D[θ[y, t], t] == 1/Pr*D[ θ[y, t], {y, 2}] - R/Pr*θ[y, t];
pde3 = D[c[y, t], t] == 1/Sc*D[c[y, t], {y, 2}]
ic = {u[y, 0] == 0, θ[y, 0] == 0, c[y, 0] == 0}
bc = {u[0, t] == 0, θ[0, t] == t, c[0, t] == t, 
     (D[u[y, t], y] /. y -> 0) == 0, (D[θ[y, t], y] /. y -> 0) == 0, 
     (D[c[y, t], y] /. y -> 0) == 0};
parm = {Gr -> 1, Gm -> 2, R -> 3, Pr -> 4, Sc -> 5};
NDSolve[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ, c}, {t, 0, 2}, {y, 0, 1}];

Update to make plot

You can make plots as follows

{usol, θsol, Csol} = 
 NDSolveValue[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ, 
   c}, {t, 0, 2}, {y, 0, 1}]
Plot3D[usol[t, y], {t, 0, 2}, {y, 0, 1}, PlotLabel -> "u solution", 
 AxesLabel -> {"time", "y", "u(t,y)"}, BaseStyle -> 14]

Plot3D[θsol[t, y], {t, 0, 2}, {y, 0, 1}, 
 PlotLabel -> "θ solution", 
 AxesLabel -> {"time", "y", "θ(t,y)"}, BaseStyle -> 14]

and so on. To plot one solution against the other, you could look into ParametricPlot and ParametricPlot3D. I am not sure what you want to plot, but there are examples how to use these commands in help.

Update to plot with varying time

like when varying t other parameters are constant.

You could use Manipulate for this. Here is an example for u(y,t) which you can the code for the other solutions. To play with varying the other parameters, you just make a slider for each one and then solve again whenever any parameter changes.

Manipulate[
 Grid[{
   {Row[{"Time = ", t0, " seconds"}]},
   {Plot[usol[y, t0], {y, 0, 1}, 
     AxesLabel -> {"y", "u(t,y)"}, BaseStyle -> 12, 
     ImageSize -> 300, PlotRange -> {Automatic, {-.1, .3}}, 
     PlotStyle -> Red, GridLines -> Automatic, 
     GridLinesStyle -> LightGray]
    }}],
 {{t0, 0, "time"}, 0, 2, .01, Appearance -> "Labeled"},
 TrackedSymbols :> {t0}
 ]

how I get 3 different line accordance to this 3 change value in R

You can make new Manipulate variable for R and solve new each time this is clicked. Like this

Make similar variable for the other variables. Here is the new code

Manipulate[
 Module[{u, y, t, c, θ, Gr, Gm, Pr, Sc, pde1, pde2, pde3},
  pde1 = 
   D[u[y, t], t] == 
    Gr*θ[y, t] + Gm*c[y, t] + D[u[y, t], {y, 2}];
  pde2 = 
   D[θ[y, t], t] == 
    1/Pr*D[θ[y, t], {y, 2}] - R/Pr*θ[y, t]; 
  pde3 = D[c[y, t], t] == 1/Sc*D[c[y, t], {y, 2}];
  ic = {u[y, 0] == 0, θ[y, 0] == 0, c[y, 0] == 0};
  bc = {u[0, t] == 0, θ[0, t] == t, 
    c[0, t] == t, (D[u[y, t], y] /. y -> 0) == 
     0, (D[θ[y, t], y] /. y -> 0) == 
     0, (D[c[y, t], y] /. y -> 0) == 0};
  parm = {Gr -> 1, Gm -> 2, Pr -> 4, Sc -> 5};
  {usol, θsol, Csol} = 
   NDSolveValue[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ, 
     c}, {y, 0, 1}, {t, 0, 2}];
  Grid[{
    {Row[{"Time = ", t0, " seconds"}]},
    {Row[{"R = ", R}]},
    {Plot[usol[y, t0], {y, 0, 1}, AxesLabel -> {"y", "u(t,y)"}, 
      BaseStyle -> 12, ImageSize -> 300, 
      PlotRange -> {Automatic, {-.1, .3}}, PlotStyle -> Red, 
      GridLines -> Automatic, GridLinesStyle -> LightGray]
     }}]
  ],
 {{t0, 0, "time"}, 0, 2, .01, Appearance -> "Labeled"},
 {R, {3, 2, 9}},
 TrackedSymbols :> {t0, R}
 ]

How to manipulate simulationsly all three plots for R .

One way is to make 3 plots and then use show. This could be improved but gives you the idea

Updated code

Manipulate[
 Module[{u, y, t, c, θ, Gr, Gm, Pr, Sc, pde1, pde2, pde3, R, 
   parm},
  pde1 = D[u[y, t], t] == Gr*θ[y, t] + Gm*c[y, t] + D[u[y, t], {y, 2}];
  pde2 = D[θ[y, t], t] == 1/Pr*D[θ[y, t], {y, 2}] - R/Pr*θ[y, t]; 
  pde3 = D[c[y, t], t] == 1/Sc*D[c[y, t], {y, 2}];
  ic = {u[y, 0] == 0, θ[y, 0] == 0, c[y, 0] == 0};
  bc = {u[0, t] == 0, θ[0, t] == t, 
    c[0, t] == t, (D[u[y, t], y] /. y -> 0) == 
     0, (D[θ[y, t], y] /. y -> 0) == 
     0, (D[c[y, t], y] /. y -> 0) == 0};
  parm = {Gr -> 1, Gm -> 2, R -> 3, Pr -> 4, Sc -> 5};
  {usol1, θsol1, Csol1} = NDSolveValue[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ, c}, {y, 0, 1}, {t, 0, 2}];
  parm = {Gr -> 1, Gm -> 2, R -> 30, Pr -> 4, Sc -> 5};
  {usol2, θsol2, Csol2} = NDSolveValue[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ, c}, {y, 0, 1}, {t, 0, 2}];
  parm = {Gr -> 1, Gm -> 2, R -> 200, Pr -> 4, Sc -> 5};
  {usol3, θsol3, Csol3} = NDSolveValue[{pde1, pde2, pde3, ic, bc} /. parm, {u, θ, c}, {y, 0, 1}, {t, 0, 2}];
  Grid[{
    {Row[{"Time = ", t0, " seconds"}]},
    {Plot[{usol1[y, t0], usol2[y, t0], usol3[y, t0]}, {y, 0, 1}, 
      AxesLabel -> {"y", "u(t,y)"}, BaseStyle -> 12, ImageSize -> 300,
       PlotRange -> {Automatic, {-.1, .3}}, 
      PlotStyle -> {Red, Blue, Green}, GridLines -> Automatic, 
      GridLinesStyle -> LightGray, 
      PlotLegends -> {"R=3", "R=30", "R=200"}]
     }}]
  ],
 {{t0, 0, "time"}, 0, 2, .01, Appearance -> "Labeled"},
 TrackedSymbols :> {t0}
 ]