clear,clc
%Solve ODE-IVP using Standard RK-4 Method
% y' = -2ty
% y(0) = 1
t0 = 0; y0 = 1;
tEnd = 5; 
h = 0.0001; N = (tEnd - t0)/h;

% Initializing variables
T = [t0:h:tEnd]';
Y = zeros(N+1,1);
Y(1) = y0;

%% Solving using Standard RK-4 Method
for i = 1:N
  k1 = myFun(T(i),Y(i));
  k2 = myFun(T(i)+h/2, Y(i)+h/2*k1);
  k3 = myFun(T(i)+h/2, Y(i)+h/2*k2);
  k4 = myFun(T(i)+h, Y(i)+h*k2);
  
  Y(i+1) = Y(i) + h/6*(k1+2*k2+2*k3+k4);
end

%% Plot and errors
plot(T,Y);
YTrue = exp(-T.^2);
ERR = abs(YTrue - Y); 
MAX_ERR = max(ERR)