I got a Runge-Kutta method here and I solve this system using it.

So here's Runge-Kutta stuff \begin{align} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + h/2, y_n + hk_1/2) \\ k_3 &= f(t_n+h, y_n - hk_1 + 2hk_2) \\\hline y_{n+1} &= y_n + h(k_1 + 4k_2 + k_3)/6 \end{align} where $h$ is step

Here's my test system \begin{align} y'_1 &= -5y_1 - 10y_2 + 14e^{-x} \\ y'_2 &= -10y_1 - 5y_2 + 14e^{-x} \end{align} with exact solution $y_1(x) = y_2(x) = e^{-x}$

UPD: The initial condition here is $y_1(0) = y_2(0) = 1$

I need to solve it on $[0;4]$.

Well, I thought I solved it right, because I checked how the exact solution and these approximate solution plots looks like (on the left, on the right I zoomed plot until saw difference) plot and zoom

Also I checked how the plot of the difference between exact solution and approximate solution depending on step (let's call it e/h) looks like.

So $e/h$ it looks like this enter image description here

But when I checked e/h^4 dependence it looked like these enter image description here

I showed it to the teacher and she said that my solution is wrong, it's not suppose to be like these! I show my code to her asked for help but she said that she doesn't understand matlab :c

Have I really done something wrong? And if yes what I've done wrong? And if not how to prove that I'm right?

Here's my code btw

Runge-Kutta method

  function [ res_y ] = RungeKutta(dim, size, grid, step, f1,f2,y1, y2)
    h = step;
    res_y(1,1) = y1;
    res_y(2,1) = y2;
    for i=1: size

       k1(1)= f1(grid(i),y1,y2);
       k1(2)= f2(grid(i),y1,y2);

       k2(1)= f1(grid(i)+h/2, y1+h*k1(1)/2, y2+h*k1(2)/2);
       k2(2)= f2(grid(i)+h/2, y1+h*k1(1)/2, y2+h*k1(2)/2);

       k3(1)= f1(grid(i)+h, y1-h*k1(1)+2*h*k2(1), y2-h*k1(2)+2*h*k2(2));
       k3(2)= f2(grid(i)+h, y1-h*k1(1)+2*h*k2(1), y2-h*k1(2)+2*h*k2(2));

       res_y(1,i+1) = y1 + h*(k1(1) + 4*k2(1) + k3(1))/6;
       res_y(2,i+1) = y1 + h*(k1(2) + 4*k2(2) + k3(2))/6;

       y1 = res_y(1,i+1); 
       y2 = res_y(2,i+1);


Main method

    a = 0; b = 4;
    h = 0.1; % step
    t = a:h:b; %grid
    n = 2; 
    m = size(t,2);
     hold on;
         plot(t, exp(-t),'b-')
         plot(t, exp(-t),'r--')
     hold off;
    y1=1; y2 = 1;
    f1_ptr = @f1;% out = -5 * y1 - 10 * y2 + (14)*exp(-x);
    f2_ptr = @f2;% out = -10 * y1 - 5 * y2 + (14)*exp(-x);
    res_y = RungeKutta(n,m-1,t,h,f1_ptr, f2_ptr,1,1);
    hold on;
    hold off;

%e/h and e/h^4 plots

 fig_a = figure;
hold on;

counter = 0;

for h=0.001:0.01:0.1
    y1=1; y2 = 1;
    t = a:h:b;
    m = size(t,2);
    counter = counter + 1; 
    result_appr = RungeKutta(n,m-1,t,h,f1_ptr, f2_ptr,y1,y2);
    result_exact = exp(-t);
   result_difference = abs(result_appr(1, :) - result_exact);
    e1(counter) = max(result_difference);
    e2(counter) = max(result_difference);
    hh = h*h*h*h;
hold off;

fig_b = figure;
set(fig_b ,'name','e/h^4','numbertitle','off')

hold on;
hold off;

f1 function

function [ out ] = f1( x, y1, y2, alpha, beta )
if nargin == 3
    alpha = 5;
    beta = 10;

out = -alpha * y1 - beta * y2 + (alpha + beta - 1)*exp(-x);


f2 function

function [ out ] = f2( x, y1, y2, alpha, beta )
if nargin == 3
    alpha = 5;
    beta = 10;
out = -beta * y1 - alpha * y2 + (alpha + beta - 1)*exp(-x);

  • $\begingroup$ Where are the initial conditions? $\endgroup$
    – Amzoti
    Nov 23, 2013 at 14:49
  • $\begingroup$ @Amzoti oops sorry, I'll update question $\endgroup$ Nov 23, 2013 at 14:57
  • $\begingroup$ @Amzoti I updated it! $\endgroup$ Nov 23, 2013 at 15:07
  • 1
    $\begingroup$ @DanilGholtsman. When you have such a problem, after a few hours, you are no more able to see anything. In this case, ask a friend to look at it. Trust me, I starting in computing sciences 53 years ago. Does this change your results ? Please reply and post. Glad to have been able to help (hoping that this is the case). $\endgroup$ Nov 23, 2013 at 16:00
  • 2
    $\begingroup$ I need to know where the test about $e/h^4$ comes from for a third order Runge-Kutta method: can your teacher give you some written reference ? $\endgroup$ Nov 25, 2013 at 17:27

1 Answer 1


The method is a classical third-order method as computed and presented by M. Wilhelm Kutta in 1901, based on the Simpson quadrature rule.

This means that the global error will scale like $h^3$ and the local truncation error like $h^4$.

Plotting the error as computed against $h$, but in a loglog plot and for a geometric sequence of step sizes, for h = 10.^[-4:0.5:-1], gives a line that has a very nice slope $3$.

enter image description here

Another way to visualize this is to plot the difference of the computed to the exact solution divided by $h^3$. This gives an error profile plot

enter image description here

that nicely and visually converges to the leading coefficient $c(x)$ in $e(x)=c(x)h^3+O(h^4)$. Note that this is the absolute error, and that the solutions fall exponentially to zero, which explains that the error falls in the end. In a plot of the relative error one could expect the curve to continue to grow away from zero.

  • $\begingroup$ Oh, after 7 years I don't even understand what it was about what I was writing in the question. As a student I knew math a bit, but years of web development washed out almost all stuff from uni completely... But nice to see that the answer is here :) $\endgroup$ Feb 17, 2021 at 23:12

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .