//third order Runge-kutta method

t_init=input("Enter initial time>");
t_final=input("Enter final time>");
h=input("Enter time step>");
N=ceil((t_final-t_init)/h); //time steps

Y=zeros(2,N+1);
T=linspace(0,2,N+1);

Y(1,1)=1;
Y(2,1)=1;

for n=2:N+1

	k1=h*hw_2f(T(1,n-1),Y(:,n-1));
	k2=h*hw_2f(T(1,n-1)+(1/3)*h,Y(:,n-1)+(1/3)*k1);
	k3=h*hw_2f(T(1,n-1)+(2/3)*h,Y(:,n-1)+(2/3)*k2);

	Y(:,n)=Y(:,n-1)+(1/4)*k1+(3/4)*k3;

end

plot2d([T',T'],[Y(1,:)',Y(2,:)'],leg="y1@y2")