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

Решение уравнения $y' = f(x,y)$, проходящее через точку $(x, y)$ должно иметь в этой точке производную $y'$, равную $f(x,y)$. То есть оно должно касаться прямой, наклоненной под углом $а= \arctan f(x,y)$ к оси $Ox$.

Часть плоскости, каждой точке которой сопоставлен отрезок прямой так, что тангенс угла наклона его к оси $Ox$ равен значению в данной точке правой части дифференциального уравнения $y' = f(x,y)$, называется полем направлений данного дифференциального уравнения.

Построение поля направлений позволяет графически находить решения дифференциальных уравнений, не находя при этом аналитическое решение.

В 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)
}

#сеть построения поля направлений
x = seq(-20,20,0.5)
y = seq(-20,20,0.5)

plot(NULL,xlim=c(-5,5),ylim=c(-3,19),ylab="", xlab="")
for(j in x)
{
  for(k in y)
  {
    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)
  }
}

Так, для уравнения $y'=y-x^2$ получится следующее: