BokehでFitzhug-Nagumo方程式をリアルタイムプロット

以下の二つの記事の合わせ技です
FitzHugh-Nagumo方程式をオイラー法で数値計算する
https://qiita.com/mozi/items/8c09d27bc373d65dfd33
bekeh+pythonでリアルタイムプロット
https://qiita.com/mozi/items/ec7e9e1181f1778a95a5

コード全文

# myplot.py
from bokeh.plotting import figure, curdoc
from bokeh.driving import linear
import numpy as np
import random
from numpy.random import rand

p = figure(plot_width=400, plot_height=400)
p.x_range.follow="end"
p.x_range.follow_interval = 40
p.x_range.range_padding=0

r1 = p.line([], [], color="firebrick", line_width=2)
r2 = p.line([], [], color="navy", line_width=2)

ds1 = r1.data_source
ds2 = r2.data_source

#initial condition
dt=0.01
a=0.2
b=0.5
c=1.0

@linear()
def update(step):
    X=0
    Y=0
    #即席なのでforで回してしまいました、、、
    for i in range(step):
        X += (c*(X-X**3-Y))*dt
        Y += (X-b*Y+a)*dt
    ds1.data['x'].append(X)
    ds1.data['y'].append(Y)
    ds2.data['x'].append(X)
    ds2.data['y'].append(Y)
    ds1.trigger('data', ds1.data, ds1.data)
    ds2.trigger('data', ds2.data, ds2.data)

curdoc().add_root(p)

# Add a periodic callback to be run every 500 milliseconds
curdoc().add_periodic_callback(update, 1)

課題

ここの部分

    X=0
    Y=0
    #即席なのでforで回してしまいました、、、
    for i in range(step):
        X += (c*(X-X**3-Y))*dt
        Y += (X-b*Y+a)*dt

callbackを利用してリアルタイムプロットしているので、後半は非常に計算量が大きくなってしまいます
forループではなく計算結果を格納した配列の値をプロットしていくほうがいいですね