Esta es la funcion que nos sirve para la propagacion de errores para dos varialbles:

function propagacion2var(f,x0,y0,ex,ey)
    syms x
    syms y;
    y1=subs(f,{x,y},{x0,y0});
    f1=diff(f,x);
    f2=diff(f,y);
    valF1=subs(f1,{x,y},{x0,y0});
    valF2=subs(f2,{x,y},{x0,y0});
    ea=abs(valF1)*ex+abs(valF2)*ey;
    er=ea/y1;
    disp('valor de y:');
    disp(y1);
    disp('error:');
    disp(ea)
    disp('er:')
    disp(er);
end


y ejecutenlo de la siguiente forma:
>> propagacion2var('x^2+4*x*y',5,3,0.2,0.1)

y el resultado sera:
valor de y:
    85

error:
32/5

er:
32/425

el resultado nos devolvera en fracciones pero si lo dividen les dara el mismo resultado