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 Runge 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.