Question: I am trying to edit the Gray-Scott 1D equation example (last example on the page) in the odeint documentation. I have the code below and it works for the Neumann boundary conditions but I want a Dirichlet boundary condition on the left at x=0 and retain Neumann on the right at x=L. I tried to reduce the dimension of y0 by 2 as that should remove 2 ODEs. But I was just repeating the same problem with 1 less interior point and getting the same answer, which makes sense since I didn't have any place where I fixed the left boundary values. import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import time from scipy.integrate import odeint def G(u, v, f, k): return f * (1 - u) - u*v**2 def H(u, v, f, k): return -(f + k) * v + u*v**2 def grayscott1d(y, t, f, k, Du, Dv, dx): """ Differential equations for the 1-D Gray-Scott equations. The ODEs are derived using the method of lines. """ # The vectors u and v are interleaved in y. We define # views of u and v by slicing y. u = y[::2] v = y[1::2] # dydt is the return value of this function. dydt = np.empty_like(y) # Just like u and v are views of the interleaved vectors # in y, dudt and dvdt are views of the interleaved output # vectors in dydt. dudt = dydt[::2] dvdt = dydt[1::2] # Compute du/dt and dv/dt. The end points and the interior points # are handled separately. dudt[0] = G(u[0], v[0], f, k) + Du * (-2.0*u[0] + 2.0*u[1]) / dx**2 dudt[1:-1] = G(u[1:-1], v[1:-1], f, k) + Du * np.diff(u,2) / dx**2 dudt[-1] = G(u[-1], v[-1], f, k) + Du * (- 2.0*u[-1] + 2.0*u[-2]) / dx**2 dvdt[0] = H(u[0], v[0], f, k) + Dv * (-2.0*v[0] + 2.0*v[1]) / dx**2 dvdt[1:-1] = H(u[1:-1], v[1:-1], f, k) + Dv * np.diff(v,2) / dx**2 dvdt[-1] = H(u[-1], v[-1], f, k) + Dv * (-2.0*v[-1] + 2.0*v[-2]) / dx**2 return dydt rng = np.random.default_rng(0) num_discretize_x = 100 #2500 gives their solution y0 = rng.standard_normal(2*num_discretize_x) y0 = np.ones(2*num_discretize_x)*3 t = np.linspace(0, 4, 100) f = 0.024 k = 0.055 Du = 0.01 Dv = 0.005 dx = 0.01 x = np.linspace(0, 1, num_discretize_x) t0 = time.time() #sola = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx)) #t1 = time.time() solb = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx), ml=2, mu=2) print('solb shape: ', solb.shape) t2 = time.time() #print(f'No banding takes {t1-t0} s, while banding takes {t2-t1} s.') u = solb[:,::2] v = solb[:,1::2] print(u.T) t_grid, x_grid = np.meshgrid(t, x) fig = plt.figure() ax = fig.add_subplot(111, projection='3d') ax.plot_surface(t_grid, x_grid, u.T) plt.xlabel('$t$') plt.ylabel('$x$') plt.show() plt.close('all') plt.plot(t, u[:,-1], label='$u(x=L)$') plt.plot(t, v[:,-1], label='$v(x=L)$') plt.xlabel('$t$') plt.ylabel('$u$ or $v$') plt.legend(loc=0) plt.show() 

Answer: You can fix the value of u and v at x = 0 and then keep the Neumann conditions at x = L: Code: import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import time from scipy.integrate import odeint def G(u, v, f, k): return f * (1 - u) - u * v**2 def H(u, v, f, k): return -(f + k) * v + u * v**2 def grayscott1d(y, t, f, k, Du, Dv, dx, u0_fixed=1.0, v0_fixed=0.5): u, v = y[::2], y[1::2] dydt = np.empty_like(y) dudt, dvdt = dydt[::2], dydt[1::2] dudt[0], dvdt[0] = 0, 0 u[0], v[0] = u0_fixed, v0_fixed dudt[1:-1] = G(u[1:-1], v[1:-1], f, k) + Du * np.diff(u, 2) / dx**2 dudt[-1] = G(u[-1], v[-1], f, k) + Du * (-2.0 * u[-1] + 2.0 * u[-2]) / dx**2 dvdt[1:-1] = H(u[1:-1], v[1:-1], f, k) + Dv * np.diff(v, 2) / dx**2 dvdt[-1] = H(u[-1], v[-1], f, k) + Dv * (-2.0 * v[-1] + 2.0 * v[-2]) / dx**2 return dydt rng = np.random.default_rng(0) num_discretize_x = 100 y0 = rng.standard_normal(2 * num_discretize_x) y0 = np.ones(2 * num_discretize_x) * 3 t = np.linspace(0, 4, 100) f = 0.024 k = 0.055 Du = 0.01 Dv = 0.0005 dx = 0.01 x = np.linspace(0, 1, num_discretize_x) t0 = time.time() solb = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx), ml=2, mu=2) print('solb shape: ', solb.shape) t2 = time.time() u = solb[:, ::2] v = solb[:, 1::2] print(u.T) t_grid, x_grid = np.meshgrid(t, x) fig = plt.figure() ax = fig.add_subplot(111, projection='3d') ax.plot_surface(t_grid, x_grid, u.T) plt.xlabel('$t$') plt.ylabel('$x$') plt.show() plt.close('all') plt.plot(t, u[:, -1], label='$u(x=L)$') plt.plot(t, v[:, -1], label='$v(x=L)$') plt.xlabel('$t$') plt.ylabel('$u$ or $v$') plt.legend(loc=0) plt.show()  Prints solb shape: (100, 200) [[3.00000000e+00 1.00000000e+00 1.00000000e+00 ... 1.00000000e+00 1.00000000e+00 1.00000000e+00] [3.00000000e+00 1.17379393e+00 8.52324657e-01 ... 3.41682696e-01 3.41131018e-01 3.40594796e-01] [3.00000000e+00 1.34513346e+00 7.43775028e-01 ... 9.71922350e-02 9.68656847e-02 9.65485574e-02] ... [3.00000000e+00 1.77119093e+00 6.99074218e-01 ... 1.18100847e-03 1.18800490e-03 1.19504096e-03] [3.00000000e+00 1.77119093e+00 6.99074218e-01 ... 1.18100847e-03 1.18800490e-03 1.19504096e-03] [3.00000000e+00 1.77119093e+00 6.99074218e-01 ... 1.18100847e-03 1.18800490e-03 1.19504096e-03]]     Comment:  @Warren Weckesser: Instead of using the parameters u0_fixed and v0_fixed and setting u[0] and v[0] in each call of grayscott1d(), you can set the corresponding initial conditions to those values (i.e. y0[:2] = [u0_fixed, v0_fixed]. Then, since dudt[0] = 0 and dvdt[0] = 0, the solution will remain constant at those values.  Code: import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import time from scipy.integrate import odeint def G(u, v, f, k): return f * (1 - u) - u * v**2 def H(u, v, f, k): return -(f + k) * v + u * v**2 def grayscott1d(y, t, f, k, Du, Dv, dx): u, v = y[::2], y[1::2] dydt = np.empty_like(y) dudt, dvdt = dydt[::2], dydt[1::2] dudt[0], dvdt[0] = 0, 0 dudt[1:-1] = G(u[1:-1], v[1:-1], f, k) + Du * np.diff(u, 2) / dx**2 dudt[-1] = G(u[-1], v[-1], f, k) + Du * (-2.0 * u[-1] + 2.0 * u[-2]) / dx**2 dvdt[1:-1] = H(u[1:-1], v[1:-1], f, k) + Dv * np.diff(v, 2) / dx**2 dvdt[-1] = H(u[-1], v[-1], f, k) + Dv * (-2.0 * v[-1] + 2.0 * v[-2]) / dx**2 return dydt rng = np.random.default_rng(0) num_discretize_x = 100 y0 = rng.standard_normal(2 * num_discretize_x) y0 = np.ones(2 * num_discretize_x) * 3 y0[0], y0[1] = 1.0, 0.5 t = np.linspace(0, 4, 100) f = 0.024 k = 0.055 Du = 0.01 Dv = 0.005 dx = 0.01 x = np.linspace(0, 1, num_discretize_x) solb = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx), ml=2, mu=2) print('solb shape: ', solb.shape) u = solb[:, ::2] v = solb[:, 1::2] print(u.T) t_grid, x_grid = np.meshgrid(t, x) fig = plt.figure() ax = fig.add_subplot(111, projection='3d') ax.plot_surface(t_grid, x_grid, u.T) plt.xlabel('$t$') plt.ylabel('$x$') plt.show() plt.close('all') plt.plot(t, u[:, -1], label='$u(x=L)$') plt.plot(t, v[:, -1], label='$v(x=L)$') plt.xlabel('$t$') plt.ylabel('$u$ or $v$') plt.legend(loc=0) plt.show()