Jako przykładu numerycznej analizy układu równań różniczkowych użyję równań Lotki-Volterry (https://pl.wikipedia.org/wiki/Równanie_Lotki-Volterry).

rownanie1=x'[t]==(a-b*y[t])*x[t]

rownanie2=y'[t]==(c *x[t]-d)*y[t]

Rozwiązanie w postaci funkcji specjalnych istnieje, ale nie jest użyteczne. Składnia DSolve (czyli funkcji rozwiązującej r. różniczkowe zwyczajne) jest następująca:
pierwszy argument - lista równań plus ewentualnych warunków początkowych
drugi argument - lista niewiadomych funkcji
trzeci argument - zmienna po której różniczkujemy (liczymy pochodne funkcji) [z definicji r. r. zwyczajnego jest tylko jedna taka zmienna (!) , często nie podana jawnie, np r. Lotki-Volterry można zapisać bez niej jako:
x' = (a-by)x
y' = (cx-d)y
]

(* DSolve[{rownanie1,rownanie2},{x,y},t] tylko dla przykładu, rozwiązania symbolicznego nie będziemy używać *)

Aby rozwiązać r. r. zwyczajne numerycznie używamy NDSolve. Funkcja ta automatycznie zanalizuje podane równania i dobierze algorytm potrzebny do ich rozwiązania. Aby rozwiązać problem numerycznie wszystkie  symbole oraz warunki początkowe muszą mieć przypisane wartości liczbowe, np:

a=1/3

b=1/2

c=4/3

d=2

Składnia NDSolve jest następująca :

argument 1 : lista równań (nie zawierająca żadnych symboli, tylko liczby !) i warunków początkowych [ {rownanie1, rownanie2, y[0] == 4, x[0] == 4/3}]

argument 2 : lista funkcji niewiadomych [  {x,y} ] lub ewentualnie [ {x[t],y[t] }

argument 3: obowiązkowy zakres zmienności zmiennej [ {t,0,100} ]

Rozwiązanie jest w postaci listy reguł transformacyjnych zawierających nowe funkcje (InterpolatingFunction) będące rozwiązaniem danego problemu:

rozwiazanie=NDSolve[{rownanie1, rownanie2, y[0]==4,x[0]==4/3},{x,y},{t,0,100}][[1]]

Aby wydobyć z niego funkcje postępujemy jak podano niżej:

ofiara=x[t]/.rozwiazanie

drapieznik=y[t]/.rozwiazanie

W ten sposób wyprodukowaliśmy dwie nowe funkcje: ofiara[t] i drapieżnik[t], które możemy od tego momentu używać praktycznie tak samo jak wbudowanych fukcji typu sinus, np. narysować wykresy:

Plot[{drapieznik,ofiara},{t,0,50},PlotRange->All]

Plot[{Tooltip[drapieznik,"drapieznik"],Tooltip[ofiara,"ofiara"]},{t,0,50},PlotRange->All]

ParametricPlot[{drapieznik,ofiara},{t,0,100},AspectRatio->1]

Punkt krytyczny to y[0] == 2/3, x[0] == 3/2, im dalej od niego, tym większe wachania populacji.