taylor series expansion function call

I have the following code for taylor series method.

function yval=taylor(f,x0,y0,xval)
h=xval-x0;
syms x y;
y=zeros(1,6);
y(1)=y0;
y(2)=f;
for i=2:5
      y(i+1)=diff(y(i),x)+y(i)*diff(y(i),y);
end
temp=h^(0:5)./factorial(0:5);
x=x0;
y=y0;
c=eval(y);
yval=sum(c.*temp);
fprintf('The value of y at given x=%f is: %f',xval,yval)
end

I am running the above program as detailed below and getting the error. suggest me how to pass the function to the program

>> f=@(x,y)x^2+y^2
f =
  function_handle with value:
    @(x,y)x^2+y^2
>> taylor(f,1,0.8,3)
Conversion to double from function_handle is not possible.
Error in taylor (line 6)
y(2)=f;

[Walter Roberson]

y(i+1)=diff(y(i),x)+y(i)*diff(y(i),y)

Why would you calculate that? taylor series is always only with respect to one variable at a time.