Plotting a complicated function containing special functions

Question

I am new to Mathematica and I'm trying to plot the following analytical function at few time steps (say t=0.1, 0.5, 1, 10) all in one plot (similar to hold on in MATLAB).

$$\begin{align*} p(x;t)=\frac1{\sigma\sqrt{2\pi t}}&\left[\exp\left(-\frac{(x-x_0-\mu t)^2}{2\sigma^2 t}\right)+\exp\left(\frac{-4x_0\mu t-(x+x_0-\mu t)^2}{2\sigma^2 t}\right)+\right.\\&\left.\frac{2\mu}{\sigma^2}\exp\left(\frac{2\mu x}{\sigma^2}\right)\left\{1-\Phi\left(\frac{x+x_0+\mu t}{\sigma \sqrt t}\right)\right\}\right] \end{align*}$$

This function is from page 224 (with some typos corrected) of the book:

Cox, D. R., and H. D. Miller. The Theory of Stochastic Processes. Vol. 134. CRC Press, 1977.

Note that the function $\Phi[(x+x_0+\mu t)/(\sigma \sqrt t)]$ is the integral from $-\infty$ to $x$ as:

$$\Phi(x)=\frac1{\sqrt{2\pi}}\int_{-\infty}^x e^{-\frac12 y^2}\mathrm dy$$

I am aware of Mathematica's error function Erf and Normal Distribution but they integrate from $0$ to $x$ and my function is from $-\infty$ to $x$.

Here is what I have for a code as an initial one-time trial with some parameters.

ClearAll[μ, σ, t, z, z0, f1, f2, f3, f4, f5, f6, f];      
z0 = 0.01; μ = 0.1; σ = 1; t = 10;
f1 := 1/(σ*Sqrt[2*π*t])
f2[z_] := Exp[-((z - z0 - μ*t)^2)/(2*t*σ^2)];
f3[z_] := Exp[(-4*z0*μ*t - (z + z0 - μ*t)^2)/(2*t*σ^2)];
f4[z_] := (2*μ)/σ^2*Exp[(2*μ*z)/σ^2];
f5[z] := Integrate[Exp[-1/2 ((u + z0 + μ*t)/(σ Sqrt[t]))^2]/Sqrt[2 π], {u, -∞, z}]
(*f6[z_]:=Simplify[f5,t>0];*)
f = f1[z]*(f2[z] + f3[z] + f4[z]*(1 - f5[z]))

(1/(10 Sqrt[2 π]))[z] 
  (E^(1/200 (-0.4 - (-9.99 + z)^2)) + E^(- (1/200) (-10.01 + z)^2) + 
    0.2 E^(0.2 z) 
      (1 - 
         ((5. - 1.32904*10^-15 I) + 
           (0.353553 (10.01 + 1. z) Erf[0.0707107 Sqrt[(10.01 + 1. z)^2]]) / 
              Sqrt[(0.707814 + 0.0707107 z)^2])[z]))

Plot[{f[z]}, {z, z0, 10}]

which is not working (I use z for x here). The function should be real, but my answer has an imaginary part.

Answer 1

I fixed some of the function definitions and changed enforced evaluation of Integrate at the time of definition with Evaluate (otherwise the Plot takes forever as the symbolic integration has to be done over and over again).

ClearAll[μ, σ, t, z, z0, f1, f2, f3, f4, f5, f6, f];
z0 = 0.01; μ = 0.1; σ = 1; t = 10;
f1[z_] := 1/(σ*Sqrt[2*π*t])
f2[z_] := Exp[-((z - z0 - μ*t)^2)/(2*t*σ^2)];
f3[z_] := Exp[(-4*z0*μ*t - (z + z0 - μ*t)^2)/(2*t*σ^2)];
f4[z_] := (2*μ)/σ^2*Exp[(2*μ*z)/σ^2];
f5[z_] := Evaluate[Integrate[
   Exp[-1/2 ((u + z0 + μ*t)/(σ Sqrt[t]))^2]/
    Sqrt[2 π], {u, -∞, z}]]
(*f6[z_]:=Simplify[f5,t>0];*)

f[z_] := f1[z]*(f2[z] + f3[z] + f4[z]*(1 - f5[z]))

Plot[{f[z]}, {z, z0, 10}]

At least, this returns plot.