
function z=f(t,y)
  z = exp(t)+3*y
endfunction


function y=rungekutta4(t,y0)
  n=prod(size(t))
  y=zeros(n,1)
  y(1)=y0
  for i=1:n-1
    hi = t(i+1)-t(i) 
    ti = t(i)
    yi1 = y(i)
    yi2 = yi1 + f(ti, yi1) * hi/2
    yi3 = yi1 + f(ti+hi/2, yi2) * hi/2
    yi4 = yi1 + f(ti+hi/2, yi3) * hi/2
    y(i+1) = yi1 + hi * (f(ti, yi1)/6 + f(ti+hi/2, yi2)/3 + f(ti+hi/2, yi3)/3 + f(ti+hi, yi4)/6)
  end
endfunction

rungekutta4([0:1/8:1],12)

for p=3:10
  tps=[0:0.5^p:1];
  plot(tps,rungekutta4(tps,12));
end
