y0 = 1;
z0 = - 1.817942;
a = 0;
b = 100;
f1 = @(x,y,z) z;
f2 = @(x,y,z) 6*y^2 - x;
x(1) = [0];
y(1) = [y0];
z(1) = [z0];
h = 0.001;
n = round((b-a)/h);
for i=2:n;
    x(i) = a + i*h;
    k1y = f1(x(i-1),y(i-1),z(i-1));
    k1z = f2(x(i-1),y(i-1),z(i-1));
    k2y = f1(x(i-1) + h/2,y(i-1) + h/2*k1y,z(i-1) + h/2*k1z);
    k2z = f2(x(i-1) + h/2,y(i-1) + h/2*k1y,z(i-1) + h/2*k1z);
    k3y = f1(x(i-1) + h/2,y(i-1) + h/2*k2y,z(i-1) + h/2*k2z);
    k3z = f2(x(i-1) + h/2,y(i-1) + h/2*k2y,z(i-1) + h/2*k2z);
    k4y = f1(x(i-1) + h/2,y(i-1) + h/2*k3y,z(i-1) + h/2*k3z);
    k4z = f2(x(i-1) + h/2,y(i-1) + h/2*k3y,z(i-1) + h/2*k3z);
    y(i) = y(i-1) + h/6*(k1y+2*k2y+2*k3y+k4y);
    z(i) = z(i-1) + h/6*(k1z+2*k2z+2*k3z+k4z);
end
plot(x,y)
hold on