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


function y=euler(t,y0)
  n=prod(size(t))
  y=zeros(n,1)
  y(1)=y0
  for i=1:n-1
    y(i+1)=y(i)+(t(i+1)-t(i))*f(t(i),y(i))
  end
endfunction

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

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