dy/dx=f(x,y)=3x+3yでy(0)=1とします。区間[0,2]を20等分した点におけるy=y(x)の値をホイン法で求めるときホイン法のプログラムは

function [x,y]=Heun(a,b,n,y_0)
h=(b-a)/n;
x=a:n:b;
y=0*x;
y(1)=y_0;
for k=1:n
s_1=h*fval([x(k),y(k)]);
s_2=h*fval([x(k)+h,y(k)+s_1]);
y(k+1)=y(k)+(s_1+s_2)/2;
end
plot(x,y,'-g');

になるんですがMATLABではエラーが出ます
s_1=h*fval([x(k),y(k)]);この部分です

[x,y]=Heun(0,2,20,1)とうつと
??? Attempted to access x(2); index out of bounds because numel(x)=1.

エラー ==> Heun at 7
s_1=h*fval([x(k),y(k)]);と表示されます。

関数は
function z=fval(w)
x=w(1);
y=w(2);
z=3*x+3*y;としています。
何が違うのかわからないです。お願いします。