Find the normal vector at a specific point on a 3D surface defined by an equation

Consider the following code:

rng('default')

[Xs,Ys,Zs] = sphere;

X = (0.99+0.01*rand(size(Xs))).*Xs;
Y = (0.99+0.01*rand(size(Ys))).*Ys;
Z = (0.99+0.01*rand(size(Zs))).*Zs;

noVertcs = 5;

xdata = X(1:noVertcs,1:noVertcs);
ydata = Y(1:noVertcs,1:noVertcs);
zdata = Z(1:noVertcs,1:noVertcs);

ft = fittype(...
    'a*x^3 + b*y^3 + c' ... 
    ,'independent', {'x','y'} ... 
    ,'dependent',{'z'} ...
    ,'coefficients',{'a','b','c'}...
    );

[XOut, YOut, ZOut] = prepareSurfaceData(xdata, ydata, zdata);

SurFit = fit( [XOut, YOut] , ZOut , ft  );

plot(SurFit)

I need to make a code that takes x, y, z coordinates of some certain point on the surface created by SurFit, then it returns the components and the direction cosines of the normal vector at this specific point.