% 4th order Runge-Kutta Demo

function d = fcn1(x,y)
  d = y - sin(x) + cos(x);
end


h        = 0.1; % 'Enter step size;
printint = 2; % 'Enter print interval;


x      = 0;    % initial conditions
xmax   = 5;
yRK    = 1;
yeuler = yRK;

write("x ","yEuler ","yRK4 ");


while x < xmax

   yeuler = yeuler + fcn1(x, yeuler) * h;    % Euler solution

   k1 = fcn1(x,         yRK);                % 4th order Runge-Kutta solution
   k2 = fcn1(x + h / 2, yRK + h * k1 / 2);
   k3 = fcn1(x + h / 2, yRK + h * k2 / 2);
   k4 = fcn1(x + h,     yRK + h * k3);
   yRK = yRK + (k1 + 2*k2 + 2*k3 + k4)*h / 6;
  
   x   = x + h;

   write(x," ",yeuler," ",yRK);
end