# MATLAB: 4th order Runge-Kutta code that can solve for several intial conditions

4th order runge kuttaMATLAB

I created a code that solves differential equations using 4th order runge-kutta method. This code can only take one intial condition. I want to make the code such that it can accept many initial conditions input as a vector and solve for each of them and store all the results in a matrix. Please help.
``% Initial conditionst0=0;x0=0;y0=1;dt=0.1;tz=10;t_range= t0:dt:tz;x_rk=zeros(1,length(t_range));y_rk= zeros(1,length(t_range));x_rk(1)=x0;y_rk(1)=y0;for i=1:length(t_range)-1      x_rk(i+1)=rk_method_x(x_rk(i),y_rk(i),dt);   y_rk(i+1)=rk_method_y(x_rk(i),y_rk(i),dt);end%% PLotsfigure;plot (t_range, x_rk,'b-o','MarkerSize',5);hold onplot (t_range, y_rk,'y-o','MarkerSize',5);%% Functionsfunction dxdt= f(x,y)dxdt= 5*x-3*y;end function dydt=g(x,y)dydt= -6*x+2*y;endfunction x1=rk_method_x(x0,y0,dt)k11=f(x0,y0);     k12=g(x0, y0);k21=f(x0+0.5*dt*k11,y0+0.5*dt*k12);k22=g(x0+0.5*dt*k11,y0+0.5*dt*k12);k31=f(x0+0.5*dt*k21,y0+0.5*dt*k22);k32=g(x0+0.5*dt*k21,y0+0.5*dt*k22);k41=f(x0+dt*k31, y0+dt*k32);k42=g(x0+dt*k31,y0+dt*k32);x1=x0+dt*((k11/6)+((k21+k31)/3)+(k41/6));endfunction y1=rk_method_y(x0,y0,dt)k11=f(x0,y0);     k12=g(x0, y0);k21=f(x0+0.5*dt*k11,y0+0.5*dt*k12);k22=g(x0+0.5*dt*k11,y0+0.5*dt*k12);k31=f(x0+0.5*dt*k21,y0+0.5*dt*k22);k32=g(x0+0.5*dt*k21,y0+0.5*dt*k22);k41=f(x0+dt*k31, y0+dt*k32);k42=g(x0+dt*k31,y0+dt*k32);y1=y0+dt*((k12/6)+((k22+k32)/3)+(k42/6));end``

``% Initial conditionst0=0;x0=0:0.2:1;y0=1:0.2:2;dt=0.1;tz=1;t_range= t0:dt:tz;X = zeros(numel(x0),numel(t_range));Y = zeros(numel(y0),numel(t_range));x_rk=zeros(1,numel(t_range));y_rk= zeros(1,numel(t_range));for n = 1:numel(x0)    x_rk(1)=x0(n);    y_rk(1)=y0(n);    for i=1:length(t_range)-1        [x_rk(i+1), y_rk(i+1)]=rk_method(x_rk(i),y_rk(i),dt);    end    X(n,:) = x_rk;    Y(n,:) = y_rk;    %% PLots    subplot(numel(x0)/2,2,n)    plot (t_range, x_rk,'b-o',t_range, y_rk,'y-o','MarkerSize',5);    grid    title([num2str(x0(n)),'   ',num2str(y0(n))])    legend('x','y')enddisp('X = '), disp(X)disp('Y = '), disp(Y)%% Functionsfunction [dxdt, dydt]= f(x,y)dxdt =  5*x-3*y;dydt = -6*x+2*y;end function [x1, y1]=rk_method(x0,y0,dt)[k11, k12]=f(x0,y0);     [k21, k22]=f(x0+0.5*dt*k11,y0+0.5*dt*k12); [k31, k32]=f(x0+0.5*dt*k21,y0+0.5*dt*k22); [k41, k42]=f(x0+dt*k31, y0+dt*k32); x1=x0+dt*((k11/6)+((k21+k31)/3)+(k41/6));y1=y0+dt*((k12/6)+((k22+k32)/3)+(k42/6));end``