function [x, t] = RungeKuttaEx(f, x0, I, tau, b, A)

b = b(:);
d = length(x0);
s = length(b);

t = I(1):tau:I(2);
N = length(t);

x = zeros(d, N);
x(:,1) = x0;


for j=2:N
    k = zeros(d,s);
    for i=1:s
        k(:,i) = f(x(:,j-1) + tau*k*(A(i,:)'));
    end
    x(:,j) = x(:,j-1) + tau*k*b;
end
