My goal is to replicate a plot similar to this one:

Although this is a physics problem. I am basically plotting and `erf()`

function. Here is how my code is currently set up right now.

```
using SpecialFunctions, Plots
T = 10^6 # temperature of solar corona [K]
η = 10^9*T^(-3/2) # m^2 s^-1
x = range(-1e7, 1e7, step=1)
x_scaled = @. x/(2*sqrt(η))
t = [0, 1, 2, 4]
B_z(x, t) = @. erf(x/sqrt(4*η*t))
plot(x_scaled, B_z(x_scaled, 0), linestyle = :dash)
xlabel!("x")
ylabel!("B_z")
xlims!(-1e7, 1e7)
```

However, I am not quite getting the result I am looking for:

What am I missing here?