Задача Джулии - решатель Рунге-Кутты модели ФитцХью – Нагумо в частных производных

Я новичок в языке программирования Julia, поэтому не знаю, как оптимизировать код. Я слышал, что Julia должна быть быстрее по сравнению с Python, но я написал простой код Julia для решения задачи ФитцХью– Модель Нагумо, и, похоже, она не быстрее Python.

Уравнения модели ФитцХью – Нагумо:

function FHN_equation(u,v,a0,a1,d,eps,dx)
  u_t = u - u.^3 - v + laplacian(u,dx)
  v_t = eps.*(u - a1 * v - a0) + d*laplacian(v,dx)
  return u_t, v_t
end

где u и v - переменные, которые представляют собой двумерные поля (то есть двумерные массивы), а a0,a1,d,eps - параметры модели. И параметры, и переменные относятся к типу Float. dx - это параметр, который управляет разделением между точками сетки для использования функции Лапласа, которая является реализацией конечных разностей с периодическими граничными условиями.

Если один из вас, опытных программистов, Джулия, подскажет, как улучшить работу Джулии, я буду рад услышать.

Пошаговая функция Рунге-Кутте:

function uv_rk4_step(Vs,Ps, dt)
  u = Vs.u
  v = Vs.v
  a0=Ps.a0
  a1=Ps.a1
  d=Ps.d
  eps=Ps.eps
  dx=Ps.dx
  du_k1, dv_k1 =    FHN_equation(u,v,a0,a1,d,eps,dx)
  u_k1 = dt*du_k1י
  v_k1 = dt*dv_k1
  du_k2, dv_k2 =    FHN_equation((u+(1/2)*u_k1),(v+(1/2)*v_k1),a0,a1,d,eps,dx)
  u_k2 = dt*du_k2
  v_k2 = dt*dv_k2
  du_k3, dv_k3 =    FHN_equation((u+(1/2)*u_k2),(v+(1/2)*v_k2),a0,a1,d,eps,dx)
  u_k3 = dt*du_k3
  v_k3 = dt*dv_k3
  du_k4, dv_k4 =    FHN_equation((u+u_k3),(v+v_k3),a0,a1,d,eps,dx)
  u_k4 = dt*du_k4
  v_k4 = dt*dv_k4
  u_next    =   u+(1/6)*u_k1+(1/3)*u_k2+(1/3)*u_k3+(1/6)*u_k4
  v_next    =   v+(1/6)*v_k1+(1/3)*v_k2+(1/3)*v_k3+(1/6)*v_k4
  return u_next, v_next
end

И я использовал imshow () из пакета PyPlot для построения поля u.