ex[k] = ex[k] + 0.5*(hy[k-1] - hy[k]);

преобразуем в два вызова blas axpy (y <- ax + y)

ex[1:end] =  0.5*hy[0:end-1] + ex[1:end]
ex[1:end] = -0.5*hy[1:end]   + ex[1:end]