function runge_kutta_expl
T = 10; s = 2; tau = 2^(-s); K = floor(T/tau);
y = zeros(K+1,1); y(1) = 2;
for k = 1:K
    y(k+1) = y(k)+tau*Phi((k-1)*tau,y(k),tau);
end
plot(tau*(0:K),y(1:K+1),'b-o'); hold on;

function val = Phi(t,y,tau)
m = 2; alpha = [0,1/2]; beta = [0,0;1/2,0]; gamma = [0,1];
eta = zeros(m,1);
val = 0;
for ell = 1:m
    dy = 0;
    for j = 1:ell-1
        dy = dy+beta(ell,j)*eta(j);
    end   
    eta(ell) = f(t+tau*alpha(ell),y+tau*dy);    
    val = val+gamma(ell)*eta(ell);
end

function val = f(t,y)
val = -2*y+5*cos(t);
