function zd = f(t,z)
g = 9.81; c = 0.05;                    % constants for gravity, friction

x = z(1);  y = z(2);  xd = z(3);  yd = z(4);

u  = cos(x) + sin(y);
ux = -sin(x);
uy =  ...
uxx = ...
uxy = ...
uyy = ...

% write problem as 
% A [ x'' ] = [ G1 ]  ,  let   w = [ x'' ]
%   [ y'' ]   [ G2 ]               [ y'' ]

A = [ 1 + ux^2 , ... ; ... , ... ];

G1 = - g*ux - ... ;
G2 = - g*uy - ... ;

w = A\[G1;G2];    

zd = [xd ; ... ; ... ; ... ];
