function [x, y] = RungekuttaS(f, c, y0, a, b, n)

h = (b-a)/n;
x = a:h:b;
m =length(y0);
y = zeros(m, n);


y(:,1)= y0; %Anfangswert Y(:,1) an der Stelle x0


for k= 1:n
    z1 = f(x(k), y(:,k), c, n);
    z2 = f(x(k) + (h/2), y(:,k) + (h/2)* z1, c, n); 
    z3 = f(x(k) + (h/2), y(:,k) + (h/2)* z2, c, n);
    z4 = f(x(k) + h, y(:,k) + h* z3, c, n);
   
    y(:,k+1) = y(:,k) +(h/6)*(z1 + 2*z2 + 2*z3 +z4);
    x(k + 1) = x(k) +h;
   
end