(*by Nasser M. Abbasi*)
u = Manipulate[Plot3D[f[x, t, x0, c, d], {x, minX, maxX}, {t, 0.01, simTime}, 
    PlotRange -> {Automatic, {0, maxTime}, {0, maxF}}, 
                 AxesLabel -> {"x", "t", "f(x,t)"}, 
    PlotLabel -> Style["Solution to df/dt =c d(xf)/dx + D d^f/dx^2 \n Ornstein-Ehrenfest process", Small], 
    ImagePadding -> 55, ImageSize -> 500, 
    PlotPoints -> ControlActive[{30, Round[simTime/10] + 5}, 
      {30, Round[simTime/10] + 5}], Mesh -> ControlActive[10, 10], 
      MaxRecursion -> 1], 
   {{c, 0.001, "c="}, 0.001, 0.5, 0.01, AppearanceElements -> All, 
     ContinuousAction -> False}, 
   {{d, 0.5, "D="}, 0.1, 5, 0.01, AppearanceElements -> All, 
     ContinuousAction -> True}, 
   {{simTime, 0.2, "simulation time="}, 0.1, maxTime, 0.1, 
     AppearanceElements -> All, SynchronousUpdating -> True}, 
   Delimiter, {{x0, 0, "x0="}, -10, 10, 0.1, AppearanceElements -> All, 
     ContinuousAction -> True}, 
   {{maxF, 0.1, "max F="}, 0.01, 3, 0.01, AppearanceElements -> All, 
     ContinuousAction -> True}, 
   AutorunSequencing -> {{3, 80}}, 
   Initialization :> {
      minX = -30; maxX = 30; maxTime = 500; 

      f[x_, t_, x0_, c_, Diffusion_] := Module[{\[Mu], \[Sigma]}, 
         \[Mu] = x0*Exp[(-c)*t]; 
         \[Sigma] = Sqrt[(Diffusion/c)*(1 - Exp[-2*c*t])]; 
         (1/(\[Sigma]*Sqrt[2*Pi]))*Exp[(-(1/2))*((x - \[Mu])/\[Sigma])^2]
      ]
    }
]