function rk(A,x,h,y0)
i=1;
y(i,:)=y0;
m=length(A);
b=x(length(x));
while x(i)<b
a=[x(i),y(i,:)];
for l=1:m
k1(l)=eval(A{l},a);
end
a=[x(i)+h/2,y(i,:)+h/2*k1];
for l=1:l
k2(l)=eval(A{l},a);
end
a=[x(i)+h/2,y(i,:)+h/2*k2];
for l=1:m
k3(l)=eval(A{l},a);
end
a=[x(i)+h,y(i,:)+h*k3];
for l=1:m
k4(l)=eval(A{l},a);
end
y(i+1,:)=y(i,:)+h/6*(k1+2*k2+2*k3+k4);
i=i+1;
end
y,k1,k2,k3,k4