rk4demo.m

% 4th order Runge-Kutta Demo
%
%  diff equation:  dy/dx = y - SIN(x) + COS(x)        subject to y(0) = 1
%
%      file fcn1.m        function dyxy = fcn1(x,y)
%                         dydx = y - sin(x) + cos(x);
%
%  exact solution is, y(x) = SIN(x) + EXP(x)

h        = input('Enter step size:      ');
printint = input('Enter print interval: ');

fprintf('\n\n  Step size: %g   print interval: %g\n\n', h, printint);

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

fprintf('   x       yexact      y Euler   EulerError      y RK4    RK4Error\n');
fprintf('%6.2f %10.4f\n',x,yRK);

i = 0;
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;

  yexact     = sin(x) + exp(x);
  Eulererror = (yeuler - yexact) / yexact * 100;
  RKerror    = (yRK - yexact) / yexact * 100;

  i = i + 1;
  if (rem(i,printint) == 0)   % print every printint step
    fprintf('%8.3f %10.4f %12.4f %10.3f %12.4f %12.3f\n',...
             x,yexact,yeuler,Eulererror,yRK,RKerror);
    i = 0;
  end;
end;

% Or try the one step Matlab solution...

[xODE45,yODE45] = ode45('fcn1',0,5,1);

fprintf('\n        x       ode45 solution    error\n');
for i = 1:length(xODE45)
  x = xODE45(i);
  yexact = sin(x) + exp(x);
  ODE45error = (yODE45(i) - yexact) / yexact * 100;
  fprintf('%12.4f  %12.4f  %12.5f\n',x,yODE45(i),ODE45error);
end;