%{
  Solujtion of y' = f(t,y,parameters) by standard RK4
%}
clear all
clc
close all

%% INPUT
T  = 2;
h  = 0.01;
p  = {};
f  = @(t,y,p) -y + 1 + t;
y0 = 0;
t0 = 0;
hs = 1.0;


%% CALCULCULATIONS
t  = t0;
y  = y0;
ys = [t,y];
while 1
    % Updating both "t" and "y" internally
  [t,y]=RK4_1Step(t,y,f,h,p);
  
  if t>=(ys(end,1)+hs)
      % Dynamic re-allocation... very useful 
      % but very inefficient
      ys = [ys;t,y];
  end

  % Very rough stopping condition
  if t>T
      break
  end
end

figure
plot(ys(:,1),ys(:,2),'b-');
grid on