diff --git a/.gitignore b/.gitignore
index c13115217..d892b7818 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,6 +4,5 @@
 deps/deps.jl
 Manifest.toml
 *.png
-docs/build
-.vscode
-.DS_store
\ No newline at end of file
+*.gif
+docs/build
\ No newline at end of file
diff --git a/docs/src/symbolic_tutorials/mol_heat.md b/docs/src/symbolic_tutorials/mol_heat.md
index 462803051..98c6da6db 100644
--- a/docs/src/symbolic_tutorials/mol_heat.md
+++ b/docs/src/symbolic_tutorials/mol_heat.md
@@ -165,3 +165,42 @@ end
 display(plt)
 savefig("plot.png")
 ```
+
+### Adding parameters
+
+We can also build up more complicated systems with multiple dependent variables and parameters as follows
+
+```julia
+@parameters t x
+@parameters Dn, Dp
+@variables u(..) v(..)
+Dt = Differential(t)
+Dx = Differential(x)
+Dxx = Differential(x)^2
+
+eqs  = [Dt(u(t,x)) ~ Dn * Dxx(u(t,x)) + u(t,x)*v(t,x), 
+        Dt(v(t,x)) ~ Dp * Dxx(v(t,x)) - u(t,x)*v(t,x)]
+bcs = [u(0,x) ~ sin(pi*x/2),
+       v(0,x) ~ sin(pi*x/2),
+       u(t,0) ~ 0.0, Dx(u(t,1)) ~ 0.0,
+       v(t,0) ~ 0.0, Dx(v(t,1)) ~ 0.0]
+
+domains = [t ∈ IntervalDomain(0.0,1.0),
+           x ∈ IntervalDomain(0.0,1.0)]
+
+pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)],[Dn=>0.5, Dp=>2])
+discretization = MOLFiniteDifference([x=>0.1],t)
+prob = discretize(pdesys,discretization) # This gives an ODEProblem since it's time-dependent
+sol = solve(prob,Tsit5())
+
+x = (0:0.1:1)[2:end-1]
+t = sol.t
+
+using Plots
+anim = @animate for i in 1:length(t)
+       p1 = plot(x,sol.u[i][1:9],label="u, t=$(t[i])";legend=false,xlabel="x",ylabel="u",ylim=[0,1])
+       p2 = plot(x,sol.u[i][10:end],label="v, t=$(t[i])";legend=false,xlabel="x",ylabel="v",ylim=[0,1])
+       plot(p1,p2)
+end
+gif(anim, "plot.gif",fps=30)
+```
\ No newline at end of file
diff --git a/test/MOL/MOL_1D_Diffusion.jl b/test/MOL/MOL_1D_Diffusion.jl
index 6d7a04f2b..44a7c12f1 100644
--- a/test/MOL/MOL_1D_Diffusion.jl
+++ b/test/MOL/MOL_1D_Diffusion.jl
@@ -382,3 +382,29 @@ end
         @test v_exact(x_sol, t_sol[i]) ≈ sol.u[i][l-1:end] atol=0.01
     end
 end
+
+@testset "Test 11: linear diffusion, two variables, mixed BCs, with parameters" begin
+    @parameters t x
+    @parameters Dn, Dp
+    @variables u(..) v(..)
+    Dt = Differential(t)
+    Dx = Differential(x)
+    Dxx = Differential(x)^2
+
+    eqs  = [Dt(u(t,x)) ~ Dn * Dxx(u(t,x)) + u(t,x)*v(t,x), 
+            Dt(v(t,x)) ~ Dp * Dxx(v(t,x)) - u(t,x)*v(t,x)]
+    bcs = [u(0,x) ~ sin(pi*x/2),
+        v(0,x) ~ sin(pi*x/2),
+        u(t,0) ~ 0.0, Dx(u(t,1)) ~ 0.0,
+        v(t,0) ~ 0.0, Dx(v(t,1)) ~ 0.0]
+
+    domains = [t ∈ IntervalDomain(0.0,1.0),
+            x ∈ IntervalDomain(0.0,1.0)]
+
+    pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)],[Dn=>0.5, Dp=>2])
+    discretization = MOLFiniteDifference([x=>0.1],t)
+    prob = discretize(pdesys,discretization)
+    @test prob.p == [0.5,2]
+    # Make sure it can be solved
+    sol = solve(prob,Tsit5())
+end
\ No newline at end of file