title: Implementing Runge-Kutta for Systems of DE
body: <p>I'm trying to implement the <a href="http://imageshack.us/photo/my-images/189/rkforde.png/" rel="nofollow">Runge-Kutta Method for Systems of DEs</a> in MATLAB. I'm not getting the <a href="http://imageshack.us/photo/my-images/802/591hw.png/" rel="nofollow">correct answers</a>, I'm not sure if there is something wrong in the code or the commands I use to run it.</p>
  
  <p>Here's my code:</p>
  
  <pre><code>function RKSystems(a, b, m, N, alpha, f)
      h = (b - a)/N;
      t = a;
      w = zeros(1, m);
  
      for j = 1:m
          w(j) = alpha(j);
      end
  
      fprintf('t = %.2f;', t);
      for i = 1:m
          fprintf(' w(%d) = %.10f;', i, w(i));
      end
      fprintf('\n');
  
      k = zeros(4, m);
      for i = 1:N
          for j = 1:m
             k(1, j) = h*f{j}(t, w);
          end
  
          for j = 1:m
              k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1));
          end
  
          for j = 1:m
              k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2));
          end
  
          for j = 1:m
              k(4, j) = h*f{j}(t + h, w + k(3));
          end
  
          for j = 1:m
              w(j) = w(j) + (k(1, j) + 2*k(2, j) + 2*k(3, j) + k(4, j))/6;
          end
  
          t = a + i*h;
  
          fprintf('t = %.2f;', t);
          for k = 1:m
              fprintf(' w(%d) = %.10f;', k, w(k));
          end
          fprintf('\n');
  
      end 
  end
  </code></pre>
  
  <p>I'm trying to test it out on <a href="http://imageshack.us/photo/my-images/534/591hwprob.png/" rel="nofollow">this problem</a>.
  Here are my commands and output:</p>
  
  <blockquote>
    <p>>> U1 = @(t, u) 3*u(1) + 2*u(2) - (2*t^2 + 1)*exp(2*t);</p>
    
    <p>>> U2 = @(t, u) 4*u(1) + u(2) + (t^2 + 2*t - 4)*exp(2*t);</p>
    
    <p>>> a = 0; b = 1; alpha = [1 1]; m = 2; h = 0.2; N = (b - a)/h;</p>
    
    <p>>> RKSystems(a, b, m, N, alpha, {U1 U2});</p>
    
    <p>t = 0.00; w(1) = 1.0000000000; w(2) = 1.0000000000;</p>
    
    <p>t = 0.20; w(1) = 2.2930309680; w(2) = 1.6186020410;</p>
    
    <p>t = 0.40; w(1) = 5.0379629523; w(2) = 3.7300162165;</p>
    
    <p>t = 0.60; w(1) = 11.4076339762; w(2) = 9.7009491301;</p>
    
    <p>t = 0.80; w(1) = 27.0898576892; w(2) = 25.6603894354;</p>
    
    <p>t = 1.00; w(1) = 67.1832886708; w(2) = 67.6103125539;</p>
  </blockquote>