Построение изоклин для обыкновенных дифференциальных уравнений в R

Задача построения интегральных кривых дифференциальных уравнений часто решается введением изоклин. Изоклиной называется геометрическое место точек, в которых касательные к интегральным кривым имеют одно и тоже направление. Семейство изоклин дифференциального уравнения $y' = f(x,y)$ определяется уравнением
$$f(x,y)=k,$$
где $k$ — параметр. Придавая параметру $k$ близкие числовые значения, можно получить достаточно густую сеть изоклин, с помощью которых можно приближенно построить интегральные кривые дифференциального yравнения.

Чтобы построить изоклины с нанесенными на них направлениями необходимо аналитически выразить уравнения изоклин. Так, для уравнения $y'=y-x^2$, общим уравнением узоклин является $y=x^2+k$.

В R построить изоклины и направления можно следующим образом:

#функция f(x,y)
diff <- function(x,y) { return(y-x^2) }

slopefieldline <- function(x1,y1,slp,d)
{  z = slope*(d-x1)+y1
  return(z)
}

plot(NULL,xlim=c(-5,5),ylim=c(-3,20),ylab="", xlab="")

c = seq(-3,3,1) # параметры изоклин
x = seq(-20,20,0.2) # плотность отображения направлений по x
for(j in x)
{
  for(l in c)
  {
    k<-j^2+l #уравнение изоклин
    slope = diff(j,k)
    domain = seq(j-0.07,j+0.07,0.14)
    z = slopefieldline(j,k,slope,domain)
    arrows(domain[1],z[1],domain[2],z[2],length=0)
  }    
}

#изоклины
for(j in c)
  { curve(x^2+j, -10, 10,add=TRUE)}

Параметр $k$ принимает значения $-3, -2,-1, 0,1,2,3$ . В результате получается 7 изоклин, по которым мы можем построить интегральные кривые: